From 7084dcfb19b4b94ded4fffaee568ab9f809c74d9 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Fri, 7 Jun 2024 13:56:15 +0300 Subject: [PATCH 01/42] Adding flavor_transfromation from "jpkneller-version-2.0" --- python/snewpy/flavor_transformation.py | 3246 +++++++++--------------- 1 file changed, 1267 insertions(+), 1979 deletions(-) diff --git a/python/snewpy/flavor_transformation.py b/python/snewpy/flavor_transformation.py index b0f3667bb..e5e6802ab 100644 --- a/python/snewpy/flavor_transformation.py +++ b/python/snewpy/flavor_transformation.py @@ -1,1979 +1,1267 @@ -# -*- coding: utf-8 -*- -"""Supernova oscillation physics for flavors e, X, e-bar, X-bar. - -For measured mixing angles and latest global analysis results, visit -http://www.nu-fit.org/. -""" - -from abc import abstractmethod, ABC - -import numpy as np -from astropy import units as u -from astropy import constants as c - -from .neutrino import MassHierarchy, MixingParameters - - -class FlavorTransformation(ABC): - """Generic interface to compute neutrino and antineutrino survival probability.""" - - @abstractmethod - def prob_ee(self, t, E): - """Electron neutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - float or ndarray - Transition probability. - """ - pass - - @abstractmethod - def prob_ex(self, t, E): - """X -> e neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - float or ndarray - Transition probability. - """ - pass - - @abstractmethod - def prob_xx(self, t, E): - """Flavor X neutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - float or ndarray - Transition probability. - """ - pass - - @abstractmethod - def prob_xe(self, t, E): - """e -> X neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - float or ndarray - Transition probability. - """ - pass - - @abstractmethod - def prob_eebar(self, t, E): - """Electron antineutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - float or ndarray - Transition probability. - """ - pass - - @abstractmethod - def prob_exbar(self, t, E): - """X -> e antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - float or ndarray - Transition probability. - """ - pass - - @abstractmethod - def prob_xxbar(self, t, E): - """X -> X antineutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - float or ndarray - Transition probability. - """ - pass - - @abstractmethod - def prob_xebar(self, t, E): - """e -> X antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - float or ndarray - Transition probability. - """ - pass - - -class NoTransformation(FlavorTransformation): - """Survival probabilities for no oscillation case.""" - - def __init__(self): - pass - - def prob_ee(self, t, E): - """Electron neutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - - def prob_ex(self, t, E): - """X -> e neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_ee(t,E) - - def prob_xx(self, t, E): - """Flavor X neutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. + self.prob_ee(t,E)) / 2. - - def prob_xe(self, t, E): - """e -> X neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_ee(t,E)) / 2. - - def prob_eebar(self, t, E): - """Electron antineutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - - def prob_exbar(self, t, E): - """X -> e antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_eebar(t,E) - - def prob_xxbar(self, t, E): - """X -> X antineutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. + self.prob_eebar(t,E)) / 2. - - def prob_xebar(self, t, E): - """e -> X antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_eebar(t,E)) / 2. - - -class CompleteExchange(FlavorTransformation): - """Survival probabilities for the case when the electron flavors are completely exchanged with the x flavor.""" - - def __init__(self): - pass - - def prob_ee(self, t, E): - """Electron neutrino survival probability. - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 0. - - def prob_ex(self, t, E): - """X -> e neutrino transition probability. - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_ee(t,E) - - def prob_xx(self, t, E): - """Flavor X neutrino survival probability. - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. + self.prob_ee(t,E)) / 2. - - def prob_xe(self, t, E): - """e -> X neutrino transition probability. - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_ee(t,E)) / 2. - - def prob_eebar(self, t, E): - """Electron antineutrino survival probability. - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 0. - - def prob_exbar(self, t, E): - """X -> e antineutrino transition probability. - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_eebar(t,E) - - def prob_xxbar(self, t, E): - """X -> X antineutrino survival probability. - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. + self.prob_eebar(t,E)) / 2. - - def prob_xebar(self, t, E): - """e -> X antineutrino transition probability. - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_eebar(t,E)) / 2. - - -class AdiabaticMSW(FlavorTransformation): - """Adiabatic MSW effect.""" - - def __init__(self, mix_angles=None, mh=MassHierarchy.NORMAL): - """Initialize transformation matrix. - - Parameters - ---------- - mix_angles : tuple or None - If not None, override default mixing angles using tuple (theta12, theta13, theta23). - mh : MassHierarchy - MassHierarchy.NORMAL or MassHierarchy.INVERTED. - """ - if type(mh) == MassHierarchy: - self.mass_order = mh - else: - raise TypeError('mh must be of type MassHierarchy') - - if mix_angles is not None: - theta12, theta13, theta23 = mix_angles - else: - pars = MixingParameters(mh) - theta12, theta13, theta23 = pars.get_mixing_angles() - - self.De1 = float((np.cos(theta12) * np.cos(theta13))**2) - self.De2 = float((np.sin(theta12) * np.cos(theta13))**2) - self.De3 = float(np.sin(theta13)**2) - - def prob_ee(self, t, E): - """Electron neutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - if self.mass_order == MassHierarchy.NORMAL: - return self.De3 - else: - return self.De2 - - def prob_ex(self, t, E): - """X -> e neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_ee(t,E) - - def prob_xx(self, t, E): - """Flavor X neutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. + self.prob_ee(t,E)) / 2. - - def prob_xe(self, t, E): - """e -> X neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_ee(t,E)) / 2. - - def prob_eebar(self, t, E): - """Electron antineutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - if self.mass_order == MassHierarchy.NORMAL: - return self.De1 - else: - return self.De3 - - def prob_exbar(self, t, E): - """X -> e antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_eebar(t,E) - - def prob_xxbar(self, t, E): - """X -> X antineutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. + self.prob_eebar(t,E)) / 2. - - def prob_xebar(self, t, E): - """e -> X antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_eebar(t,E)) / 2. - - -class NonAdiabaticMSWH(FlavorTransformation): - """Nonadiabatic MSW effect.""" - - def __init__(self, mix_angles=None, mh=MassHierarchy.NORMAL): - """Initialize transformation matrix. - - Parameters - ---------- - mix_angles : tuple or None - If not None, override default mixing angles using tuple (theta12, theta13, theta23). - mh : MassHierarchy - MassHierarchy.NORMAL or MassHierarchy.INVERTED. - """ - if type(mh) == MassHierarchy: - self.mass_order = mh - else: - raise TypeError('mh must be of type MassHierarchy') - - if mix_angles is not None: - theta12, theta13, theta23 = mix_angles - else: - pars = MixingParameters(mh) - theta12, theta13, theta23 = pars.get_mixing_angles() - - self.De1 = float((np.cos(theta12) * np.cos(theta13))**2) - self.De2 = float((np.sin(theta12) * np.cos(theta13))**2) - self.De3 = float(np.sin(theta13)**2) - - def prob_ee(self, t, E): - """Electron neutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return self.De2 - - def prob_ex(self, t, E): - """X -> e neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_ee(t,E) - - def prob_xx(self, t, E): - """Flavor X neutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. + self.prob_ee(t,E)) / 2. - - def prob_xe(self, t, E): - """e -> X neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_ee(t,E)) / 2. - - def prob_eebar(self, t, E): - """Electron antineutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return self.De1 - - def prob_exbar(self, t, E): - """X -> e antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_eebar(t,E) - - def prob_xxbar(self, t, E): - """X -> X antineutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. + self.prob_eebar(t,E)) / 2. - - def prob_xebar(self, t, E): - """e -> X antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_eebar(t,E)) / 2. - - -class TwoFlavorDecoherence(FlavorTransformation): - """Star-earth transit survival probability: two flavor case.""" - - def __init__(self, mix_angles=None, mh=MassHierarchy.NORMAL): - """Initialize transformation matrix. - - Parameters - ---------- - mix_angles : tuple or None - If not None, override default mixing angles using tuple (theta12, theta13, theta23). - mh : MassHierarchy - MassHierarchy.NORMAL or MassHierarchy.INVERTED. - """ - if type(mh) == MassHierarchy: - self.mass_order = mh - else: - raise TypeError('mh must be of type MassHierarchy') - - if mix_angles is not None: - theta12, theta13, theta23 = mix_angles - else: - pars = MixingParameters(mh) - theta12, theta13, theta23 = pars.get_mixing_angles() - - self.De1 = float((np.cos(theta12) * np.cos(theta13))**2) - self.De2 = float((np.sin(theta12) * np.cos(theta13))**2) - self.De3 = float(np.sin(theta13)**2) - - def prob_ee(self, t, E): - """Electron neutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - """ - if self.mass_order == MassHierarchy.NORMAL: - return (self.De2+self.De3)/2. - else: - return self.De2 - - def prob_ex(self, t, E): - """X -> e neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_ee(t,E) - - def prob_xx(self, t, E): - """Flavor X neutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. + self.prob_ee(t,E)) / 2. - - def prob_xe(self, t, E): - """e -> X neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_ee(t,E)) / 2. - - def prob_eebar(self, t, E): - """Electron antineutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - if self.mass_order == MassHierarchy.NORMAL: - return self.De1 - else: - return (self.De1+self.De3)/2 - - def prob_exbar(self, t, E): - """X -> e antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_eebar(t,E) - - def prob_xxbar(self, t, E): - """X -> X antineutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. + self.prob_eebar(t,E)) / 2. - - def prob_xebar(self, t, E): - """e -> X antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_eebar(t,E)) / 2. - - -class ThreeFlavorDecoherence(FlavorTransformation): - """Star-earth transit survival probability: three flavor case.""" - - def __init__(self): - pass - - def prob_ee(self, t, E): - """Electron neutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - """ - return 1./3. - - def prob_ex(self, t, E): - """X -> e neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_ee(t,E) - - def prob_xx(self, t, E): - """Flavor X neutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. + self.prob_ee(t,E)) / 2. - - def prob_xe(self, t, E): - """e -> X neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_ee(t,E)) / 2. - - def prob_eebar(self, t, E): - """Electron antineutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1./3. - - def prob_exbar(self, t, E): - """X -> e antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_eebar(t,E) - - def prob_xxbar(self, t, E): - """X -> X antineutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. + self.prob_eebar(t,E)) / 2. - - def prob_xebar(self, t, E): - """e -> X antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_eebar(t,E)) / 2. - - -class NeutrinoDecay(FlavorTransformation): - """Decay effect, where the heaviest neutrino decays to the lightest - neutrino. For a description and typical parameters, see A. de Gouvêa et al., - PRD 101:043013, 2020, arXiv:1910.01127. - """ - def __init__(self, mix_angles=None, mass=1*u.eV/c.c**2, tau=1*u.day, dist=10*u.kpc, mh=MassHierarchy.NORMAL): - """Initialize transformation matrix. - - Parameters - ---------- - mix_angles : tuple or None - If not None, override default mixing angles using tuple (theta12, theta13, theta23). - mass : astropy.units.quantity.Quantity - Mass of the heaviest neutrino; expect in eV/c^2. - tau : astropy.units.quantity.Quantity - Lifetime of the heaviest neutrino. - dist : astropy.units.quantity.Quantity - Distance to the supernova. - mh : MassHierarchy - MassHierarchy.NORMAL or MassHierarchy.INVERTED. - """ - if type(mh) == MassHierarchy: - self.mass_order = mh - else: - raise TypeError('mh must be of type MassHierarchy') - - if mix_angles is not None: - theta12, theta13, theta23 = mix_angles - else: - pars = MixingParameters(mh) - theta12, theta13, theta23 = pars.get_mixing_angles() - - self.De1 = float((np.cos(theta12) * np.cos(theta13))**2) - self.De2 = float((np.sin(theta12) * np.cos(theta13))**2) - self.De3 = float(np.sin(theta13)**2) - - self.m = mass - self.tau = tau - self.d = dist - - def gamma(self, E): - """Decay width of the heaviest neutrino mass. - - Parameters - ---------- - E : float - Energy of the nu3. - - Returns - ------- - Gamma : float - Decay width of the neutrino mass, in units of 1/length. - - :meta private: - """ - return self.m*c.c / (E*self.tau) - - def prob_ee(self, t, E): - """Electron neutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - # NMO case. - if self.mass_order == MassHierarchy.NORMAL: - pe_array = self.De1*(1-np.exp(-self.gamma(E)*self.d)) + \ - self.De3*np.exp(-self.gamma(E)*self.d) - # IMO case. - else: - pe_array = self.De2*np.exp(-self.gamma(E)*self.d) + \ - self.De3*(1-np.exp(-self.gamma(E)*self.d)) - return pe_array - - def prob_ex(self, t, E): - """X -> e neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - # NMO case. - if self.mass_order == MassHierarchy.NORMAL: - return self.De1 + self.De3 - # IMO case. - else: - return self.De1 + self.De2 - - def prob_xx(self, t, E): - """Flavor X neutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_ex(t,E) / 2. - - def prob_xe(self, t, E): - """e -> X neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_ee(t,E)) / 2. - - def prob_eebar(self, t, E): - """Electron antineutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return self.De3 - - def prob_exbar(self, t, E): - """X -> e antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - # NMO case. - if self.mass_order == MassHierarchy.NORMAL: - pxbar_array = self.De1*(1-np.exp(-self.gamma(E)*self.d)) + \ - self.De2 + self.De3*np.exp(-self.gamma(E)*self.d) - # IMO case. - else: - pxbar_array = self.De1 + self.De2*np.exp(-self.gamma(E)*self.d) + \ - self.De3*(1-np.exp(-self.gamma(E)*self.d)) - return pxbar_array - - def prob_xxbar(self, t, E): - """X -> X antineutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_exbar(t,E) / 2. - - def prob_xebar(self, t, E): - """e -> X antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_eebar(t,E)) / 2. - - -class AdiabaticMSWes(FlavorTransformation): - """A four-neutrino mixing prescription. The assumptions used are that: - - 1. the fourth neutrino mass is the heaviest but not so large that the electron-sterile resonances - are inside the neutrinosphere; - 2. the “outer” or H' electron-sterile MSW resonance is adiabatic; - 3. the “inner” or H'' electron-sterile MSW resonance (where the electron fraction = 1/3) is non-adiabatic. - - For further insight see, for example, Esmaili, Peres, and Serpico, Phys. Rev. D 90, 033013 (2014). - """ - def __init__(self, mix_angles, mh=MassHierarchy.NORMAL): - """Initialize transformation matrix. - - Parameters - ---------- - mix_angles : tuple - Values for mixing angles (theta12, theta13, theta23, theta14). - mh : MassHierarchy - MassHierarchy.NORMAL or MassHierarchy.INVERTED. - """ - if type(mh) == MassHierarchy: - self.mass_order = mh - else: - raise TypeError('mh must be of type MassHierarchy') - - theta12, theta13, theta23, theta14 = mix_angles - - self.De1 = float((np.cos(theta12) * np.cos(theta13) * np.cos(theta14))**2) - self.De2 = float((np.sin(theta12) * np.cos(theta13) * np.cos(theta14))**2) - self.De3 = float((np.sin(theta13) * np.cos(theta14))**2) - self.De4 = float((np.sin(theta14))**2) - self.Ds1 = float((np.cos(theta12) * np.cos(theta13) * np.sin(theta14))**2) - self.Ds2 = float((np.sin(theta12) * np.cos(theta13) * np.sin(theta14))**2) - self.Ds3 = float((np.sin(theta13) * np.sin(theta14))**2) - self.Ds4 = float((np.cos(theta14))**2) - - def prob_ee(self, t, E): - """e -> e neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return self.De4 - - def prob_ex(self, t, E): - """x -> e neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - if self.mass_order == MassHierarchy.NORMAL: - return self.De1 + self.De2 - else: - return self.De1 + self.De3 - - def prob_xx(self, t, E): - """x -> x neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - if self.mass_order == MassHierarchy.NORMAL: - return ( 2 - self.De1 - self.De2 - self.Ds1 - self.Ds2 ) / 2 - else: - return ( 2 - self.De1 - self.De3 - self.Ds1 - self.Ds3 ) / 2 - - def prob_xe(self, t, E): - """e -> x neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return ( 1 - self.De4 - self.Ds4 )/2 - - def prob_eebar(self, t, E): - """e -> e antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - if self.mass_order == MassHierarchy.NORMAL: - return self.De1 - else: - return self.De3 - - def prob_exbar(self, t, E): - """x -> e antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - if self.mass_order == MassHierarchy.NORMAL: - return self.De3 + self.De4 - else: - return self.De2 + self.De4 - - def prob_xxbar(self, t, E): - """x -> x antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - if self.mass_order == MassHierarchy.NORMAL: - return ( 2 - self.De3 - self.De4 - self.Ds3 - self.Ds4 ) / 2 - else: - return ( 2 - self.De2 - self.De4 - self.Ds2 - self.Ds4 ) / 2 - - def prob_xebar(self, t, E): - """e -> x antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - if self.mass_order == MassHierarchy.NORMAL: - return ( 1 - self.De1 - self.Ds1 ) / 2 - else: - return ( 1 - self.De3 - self.Ds3 ) / 2 - - -class NonAdiabaticMSWes(FlavorTransformation): - """A four-neutrino mixing prescription. The assumptions used are that: - - 1. the fourth neutrino mass is the heaviest but not so large that the electron-sterile resonances - are inside the neutrinosphere; - 2. the “outer” or H' electron-sterile MSW resonance is non-adiabatic; - 3. the “inner” or H'' electron-sterile MSW resonance (where the electron fraction = 1/3) is non-adiabatic. - - For further insight see, for example, Esmaili, Peres, and Serpico, Phys. Rev. D 90, 033013 (2014). - """ - def __init__(self, mix_angles, mh=MassHierarchy.NORMAL): - """Initialize transformation matrix. - - Parameters - ---------- - mix_angles : tuple - Values for mixing angles (theta12, theta13, theta23, theta14). - mh : MassHierarchy - MassHierarchy.NORMAL or MassHierarchy.INVERTED. - """ - if type(mh) == MassHierarchy: - self.mass_order = mh - else: - raise TypeError('mh must be of type MassHierarchy') - - theta12, theta13, theta23, theta14 = mix_angles - - self.De1 = float((np.cos(theta12) * np.cos(theta13) * np.cos(theta14))**2) - self.De2 = float((np.sin(theta12) * np.cos(theta13) * np.cos(theta14))**2) - self.De3 = float((np.sin(theta13) * np.cos(theta14))**2) - self.De4 = float((np.sin(theta14))**2) - self.Ds1 = float((np.cos(theta12) * np.cos(theta13) * np.sin(theta14))**2) - self.Ds2 = float((np.sin(theta12) * np.cos(theta13) * np.sin(theta14))**2) - self.Ds3 = float((np.sin(theta13) * np.sin(theta14))**2) - self.Ds4 = float((np.cos(theta14))**2) - - def prob_ee(self, t, E): - """e -> e neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - if self.mass_order == MassHierarchy.NORMAL: - return self.De3 - else: - return self.De2 - - def prob_ex(self, t, E): - """x -> e neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - if self.mass_order == MassHierarchy.NORMAL: - return self.De1 + self.De2 - else: - return self.De1 + self.De3 - - def prob_xx(self, t, E): - """x -> x neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - if self.mass_order == MassHierarchy.NORMAL: - return ( 2 - self.De1 - self.De2 - self.Ds1 - self.Ds2 ) / 2 - else: - return ( 2 - self.De1 - self.De3 - self.Ds1 - self.Ds3 ) / 2 - - def prob_xe(self, t, E): - """e -> x neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - if self.mass_order == MassHierarchy.NORMAL: - return (1 - self.De3 - self.Ds3)/2 - else: - return (1 - self.De2 - self.Ds2) / 2 - - def prob_eebar(self, t, E): - """e -> e antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - if self.mass_order == MassHierarchy.NORMAL: - return self.De1 - else: - return self.De3 - - def prob_exbar(self, t, E): - """x -> e antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - if self.mass_order == MassHierarchy.NORMAL: - return self.De2 + self.De3 - else: - return self.De1 + self.De2 - - def prob_xxbar(self, t, E): - """x -> x antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - if self.mass_order == MassHierarchy.NORMAL: - return ( 2 - self.De2 - self.De3 - self.Ds2 - self.Ds3 ) / 2 - else: - return ( 2 - self.De1 - self.De2 - self.Ds1 - self.Ds2 ) / 2 - - def prob_xebar(self, t, E): - """e -> x antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - if self.mass_order == MassHierarchy.NORMAL: - return ( 1 - self.De1 - self.Ds1 ) / 2 - else: - return ( 1 - self.De3 - self.Ds3 ) / 2 - - -class QuantumDecoherence(FlavorTransformation): - """Quantum Decoherence, where propagation in vacuum leads to equipartition - of states. For a description and typical parameters, see M. V. dos Santos et al., - 2023, arXiv:2306.17591. - """ - def __init__(self, mix_angles=None, Gamma3=1e-27*u.eV, Gamma8=1e-27*u.eV, dist=10*u.kpc, n=0, E0=10*u.MeV, mh=MassHierarchy.NORMAL): - """Initialize transformation matrix. - - Parameters - ---------- - mix_angles : tuple or None - If not None, override default mixing angles using tuple (theta12, theta13, theta23). - Gamma3 : astropy.units.quantity.Quantity - Quantum decoherence parameter; expect in eV. - Gamma8 : astropy.units.quantity.Quantity - Quantum decoherence parameter; expect in eV. - dist : astropy.units.quantity.Quantity - Distance to the supernova. - n : float - Exponent of power law for energy dependent quantum decoherence parameters, - i.e. Gamma = Gamma0*(E/E0)**n. If not specified, it is taken as zero. - E0 : astropy.units.quantity.Quantity - Reference energy in the power law Gamma = Gamma0*(E/E0)**n. If not specified, - it is taken as 10 MeV. Note that if n = 0, quantum decoherence parameters are independent - of E0. - mh : MassHierarchy - MassHierarchy.NORMAL or MassHierarchy.INVERTED. - """ - if type(mh) == MassHierarchy: - self.mass_order = mh - else: - raise TypeError('mh must be of type MassHierarchy') - - if mix_angles is not None: - theta12, theta13, theta23 = mix_angles - else: - pars = MixingParameters(mh) - theta12, theta13, theta23 = pars.get_mixing_angles() - - self.De1 = float((np.cos(theta12) * np.cos(theta13))**2) - self.De2 = float((np.sin(theta12) * np.cos(theta13))**2) - self.De3 = float(np.sin(theta13)**2) - - self.Gamma3 = (Gamma3 / (c.hbar.to('eV s') * c.c)).to('1/kpc') - self.Gamma8 = (Gamma8 / (c.hbar.to('eV s') * c.c)).to('1/kpc') - self.d = dist - self.n = n - self.E0 = E0 - - def P11(self, E): - """Survival probability of state nu1 in vacuum. - - Parameters - ---------- - E : float - Energy. - - Returns - ------- - P11 : float - Survival probability of state nu1 in vacuum. - - :meta private: - """ - return 1/3 + 1/2 * np.exp(-(self.Gamma3 * (E/self.E0)**self.n + self.Gamma8 * (E/self.E0)**self.n / 3) * self.d) + 1/6 * np.exp(-self.Gamma8 * (E/self.E0)**self.n * self.d) - - def P21(self, E): - """Transition probability from the state nu2 to nu1 in vacuum. - - Parameters - ---------- - E : float - Energy. - - Returns - ------- - P21 : float - Transition probability from the state nu2 to nu1 in vacuum. - Note that P21 = P12. - - :meta private: - """ - return 1/3 - 1/2 * np.exp(-(self.Gamma3 * (E/self.E0)**self.n + self.Gamma8 * (E/self.E0)**self.n / 3) * self.d) + 1/6 * np.exp(-self.Gamma8 * (E/self.E0)**self.n * self.d) - - def P22(self, E): - """Survival probability of state nu2 in vacuum. - - Parameters - ---------- - E : float - Energy. - - Returns - ------- - P21 : float - Survival probability of state nu2 in vacuum. - - :meta private: - """ - return self.P11(E) - - - def P31(self, E): - """Transition probability from the state nu3 to nu1 in vacuum. - - Parameters - ---------- - E : float - Energy. - - Returns - ------- - P31 : float - Transition probability from the state nu3 to nu1 in vacuum. - Note that P31 = P13. - - :meta private: - """ - return 1/3 - 1/3 * np.exp(-self.Gamma8 * (E/self.E0)**self.n * self.d) - - def P32(self, E): - """Transition probability from the state nu3 to nu2 in vacuum. - - Parameters - ---------- - E : float - Energy. - - Returns - ------- - P32 : float - Transition probability from the state nu3 to nu2 in vacuum. - Note that P32 = P23. - - :meta private: - """ - return self.P31(E) - - def P33(self, E): - """Survival probability of state nu3 in vacuum. - - Parameters - ---------- - E : float - Energy. - - Returns - ------- - P33 : float - Survival probability of state nu3 in vacuum. - - :meta private: - """ - return 1/3 + 2/3 * np.exp(-self.Gamma8 * (E/self.E0)**self.n * self.d) - - def prob_ee(self, t, E): - """Electron neutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - # NMO case. - if self.mass_order == MassHierarchy.NORMAL: - pe_array = self.P31(E)*self.De1 + self.P32(E)*self.De2 + self.P33(E)*self.De3 - # IMO case. - else: - pe_array = self.P22(E)*self.De2 + self.P21(E)*self.De1 + self.P32(E)*self.De3 - return pe_array - - def prob_ex(self, t, E): - """X -> e neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_ee(t,E) - - def prob_xx(self, t, E): - """Flavor X neutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_ex(t,E) / 2. - - def prob_xe(self, t, E): - """e -> X neutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_ee(t,E)) / 2. - - def prob_eebar(self, t, E): - """Electron antineutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - # NMO case. - if self.mass_order == MassHierarchy.NORMAL: - pe_array = self.P11(E)*self.De1 + self.P21(E)*self.De2 + self.P31(E)*self.De3 - # IMO case. - else: - pe_array = self.P31(E)*self.De1 + self.P32(E)*self.De2 + self.P33(E)*self.De3 - return pe_array - - def prob_exbar(self, t, E): - """X -> e antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_eebar(t,E) - - def prob_xxbar(self, t, E): - """X -> X antineutrino survival probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_exbar(t,E) / 2. - - def prob_xebar(self, t, E): - """e -> X antineutrino transition probability. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_eebar(t,E)) / 2. +# -*- coding: utf-8 -*- +"""Supernova oscillation physics + +For measured mixing angles and latest global analysis results, visit +http://www.nu-fit.org/. +""" + +from abc import abstractmethod, ABC + +import numpy as np +from astropy import units as u +from astropy import constants as c +from astropy.coordinates import AltAz + +from .neutrino import MassHierarchy, ThreeFlavor, FourFlavor + +############################################################################### + +class FlavorTransformation(ABC): + """Generic interface to compute neutrino and antineutrino survival probability.""" + + @abstractmethod + def get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities. + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + p : a N x N array, or an array of N x N arrays + where N is either 6 or 8 + """ + pass + +############################################################################### + +class ThreeFlavorTransformation(FlavorTransformation): + """Base class defining common data and methods for all three flavor transformations""" + + def __init__(self, mix_params): + """Initialize flavor transformation + + Parameters + ---------- + mix_params : ThreeFlavorMixingParameters instance or None + """ + if mix_params == None: + self.mix_params = ThreeFlavorMixingParameters(MassHierarchy.NORMAL) + else: + self.mix_params = mix_params + + + def Pmf_HighDensityLimit(self): + """ The probability that a given flavor state is a particular matter state in the + infinite density limit. + + Returns + ------- + 6 x 6 matrix + """ + PmfHDL = np.zeros((6,6)) + + NU_E, NU_MU, NU_TAU, NU_E_BAR, NU_MU_BAR, NU_TAU_BAR = \ + ThreeFlavor.NU_E, ThreeFlavor.NU_MU, ThreeFlavor.NU_TAU, \ + ThreeFlavor.NU_E_BAR, ThreeFlavor.NU_MU_BAR, ThreeFlavor.NU_TAU_BAR + + M2 = np.zeros((6,6)) + M2[1,1] = self.mix_params.dm21_2.value + M2[2,2] = self.mix_params.dm31_2.value + M2[1+3,1+3] = self.mix_params.dm21_2.value + M2[2+3,2+3] = self.mix_params.dm31_2.value + + U = self.mix_params.VacuumMixingMatrix() + + HV = U @ M2 @ np.conjugate(np.transpose(U)) + + T = np.real( ( HV[NU_MU,NU_MU] + HV[NU_TAU,NU_TAU] ) / 2 ) + D = np.abs( HV[NU_MU,NU_TAU] )**2 - np.abs( HV[NU_MU,NU_MU]-HV[NU_TAU,NU_TAU] )**2 / 4 + + Tbar = np.real( ( HV[NU_MU_BAR,NU_MU_BAR] + HV[NU_TAU_BAR,NU_TAU_BAR] ) / 2 ) + Dbar = np.abs( HV[NU_MU_BAR,NU_TAU_BAR] )**2 \ + -np.abs( HV[NU_MU_BAR,NU_MU_BAR]-HV[NU_TAU_BAR,NU_TAU_BAR] )**2 / 4 + + """ The NMO case. Matter state 3 is the electron neutrino, matter state 1bar is the electron + antineutrino. + """ + if self.mix_params.mass_order == MassHierarchy.NORMAL: + k1 = T - np.sqrt(D) + k2 = T + np.sqrt(D) + k2bar = Tbar - np.sqrt(Dbar) + k3bar = Tbar + np.sqrt(Dbar) + + PmfHDL[0,NU_MU] = ( np.real(HV[NU_TAU,NU_TAU]) - k1 ) / ( k2 - k1 ) + PmfHDL[0,NU_TAU] = ( np.real(HV[NU_MU,NU_MU]) - k1 ) / ( k2 - k1 ) + PmfHDL[1,NU_MU] = ( np.real(HV[NU_TAU,NU_TAU]) - k2 ) / ( k1 - k2 ) + PmfHDL[1,NU_TAU] = ( np.real(HV[NU_MU,NU_MU]) - k2 ) / ( k1 - k2 ) + PmfHDL[2,NU_E] = 1 + + PmfHDL[3,NU_E_BAR] = 1 + PmfHDL[4,NU_MU_BAR] = ( np.real(HV[NU_TAU_BAR,NU_TAU_BAR]) - k2bar ) / ( k3bar - k2bar ) + PmfHDL[4,NU_TAU_BAR] = ( np.real(HV[NU_MU_BAR,NU_MU_BAR]) - k2bar ) / ( k3bar - k2bar ) + PmfHDL[5,NU_MU_BAR] = ( np.real(HV[NU_TAU_BAR,NU_TAU_BAR]) - k3bar ) / ( k2bar - k3bar ) + PmfHDL[5,NU_TAU_BAR] = ( np.real(HV[NU_MU_BAR,NU_MU_BAR]) - k3bar ) / ( k2bar - k3bar ) + + """ The IMO case. Matter state 2 is the electron neutrino, matter state 3bar is the electron + antineutrino. + """ + if self.mix_params.mass_order == MassHierarchy.INVERTED: + k1 = T + np.sqrt(D) + k3 = T - np.sqrt(D) + k1bar = Tbar - np.sqrt(Dbar) + k2bar = Tbar + np.sqrt(Dbar) + + PmfHDL[0,NU_MU] = ( np.real(HV[NU_TAU,NU_TAU]) - k1 ) / ( k3 - k1 ) + PmfHDL[0,NU_TAU] = ( np.real(HV[NU_MU,NU_MU]) - k1 ) / ( k3 - k1 ) + PmfHDL[1,NU_E] = 1 + PmfHDL[2,NU_MU] = ( np.real(HV[NU_TAU,NU_TAU]) - k3 ) / ( k1 - k3 ) + PmfHDL[2,NU_TAU] = ( np.real(HV[NU_MU,NU_MU]) - k3 ) / ( k1 - k3 ) + + PmfHDL[3,NU_MU_BAR] = ( np.real(HV[NU_TAU_BAR,NU_TAU_BAR]) - k1bar ) / ( k2bar - k1bar ) + PmfHDL[3,NU_TAU_BAR] = ( np.real(HV[NU_MU_BAR,NU_MU_BAR]) - k1bar ) / ( k2bar - k1bar ) + PmfHDL[4,NU_MU_BAR] = ( np.real(HV[NU_TAU_BAR,NU_TAU_BAR]) - k2bar ) / ( k1bar - k2bar ) + PmfHDL[4,NU_TAU_BAR] = ( np.real(HV[NU_MU_BAR,NU_MU_BAR]) - k2bar ) / ( k1bar - k2bar ) + PmfHDL[5,NU_E_BAR] = 1 + + return PmfHDL + +############################################################################### + +class FourFlavorTransformation: + """Base class defining common data and method for all four flavor transformations""" + + def __init__(self, mix_params): + """Initialize flavor transformation + + Parameters + ---------- + mix_params : FourFlavorMixingParameters instance + """ + if mix_params == None: + self.mix_params = FourFlavorMixingParameters(MassHierarchy.NORMAL) + else: + self.mix_params = mix_params + + def Pmf_HighDensityLimit(self): + """ The probability that a given flavor state is 'detected' as a particular matter state in the + infinite density limit. This method assumes the 4th matter state is the + heaviest and that the electron fraction remains larger than 1/3 + + Returns + ------- + 8 x 8 matrix + """ + PmfHDL = np.zeros((8,8)) + + NU_E, NU_MU, NU_TAU, NU_S, NU_E_BAR, NU_MU_BAR, NU_TAU_BAR, NU_S_BAR = \ + FourFlavor.NU_E, FourFlavor.NU_MU, FourFlavor.NU_TAU, FourFlavor.NU_S, \ + FourFlavor.NU_E_BAR, FourFlavor.NU_MU_BAR, FourFlavor.NU_TAU_BAR, FourFlavor.NU_S_BAR + + M2 = np.zeros((8,8)) + M2[1,1] = self.mix_params.dm21_2.value + M2[2,2] = self.mix_params.dm31_2.value + M2[3,3] = self.mix_params.dm41_2.value + M2[1+4,1+4] = self.mix_params.dm21_2.value + M2[2+4,2+4] = self.mix_params.dm31_2.value + M2[3+4,3+4] = self.mix_params.dm41_2.value + + U = VacuumMixingMatrix() + + HV = np.zeros((6,6,len(E))) + for m in range(len(E)): + HV[:,:,m] = U @ M2 @ np.conjugate(np.transpose(U)) / ( 2 * E[m] ) + + T = np.real( (HV[NU_MU,NU_MU] + HV[NU_TAU,NU_TAU]) / 2 ) + D = np.abs( HV[NU_MU,NU_TAU] )**2 - np.abs( HV[NU_MU,NU_MU]-HV[NU_TAU,NU_TAU] )**2 / 4 + + Tbar = np.real( ( HV[NU_MU_BAR,NU_MU_BAR] + HV[NU_TAU_BAR,NU_TAU_BAR] ) / 2 ) + Dbar = np.abs( HV[NU_MU_BAR,NU_TAU_BAR] )**2 \ + -np.abs( HV[NU_MU_BAR,NU_MU_BAR]-HV[NU_TAU_BAR,NU_TAU_BAR] )**2 / 4 + + U = Vacuum_Mixing_Matrix() + + """ The NMO case. Matter state 4 is the electron neutrino, matter state 1bar is the electron + antineutrino. Matter state 3 is the sterile neutrino, matter state 2bar is the sterile + antineutrino. + """ + if self.mix_params.mass_order == MassHierarchy.NORMAL: + k1 = T - np.sqrt(D) + k2 = T + np.sqrt(D) + k3bar = Tbar - np.sqrt(Dbar) + k4bar = Tbar + np.sqrt(Dbar) + + PmfHDL[0,NU_MU] = ( np.real(HV[NU_TAU,NU_TAU]) - k1 ) / ( k2 - k1 ) + PmfHDL[0,NU_TAU] = ( np.real(HV[NU_MU,NU_MU]) - k1 ) / ( k2 - k1 ) + PmfHDL[1,NU_MU] = ( np.real(HV[NU_TAU,NU_TAU]) - k2 ) / ( k1 - k2 ) + PmfHDL[1,NU_TAU] = ( np.real(HV[NU_MU,NU_MU]) - k2 ) / ( k1 - k2 ) + PmfHDL[2,NU_S] = 1 + PmfHDL[3,NU_E] = 1 + + PmfHDL[4,NU_E_BAR] = 1 + PmfHDL[5,NU_S_BAR] = 1 + PmfHDL[6,NU_MU_BAR] = ( np.real(HV[NU_TAU_BAR,NU_TAU_BAR]) - k3bar ) / ( k4bar - k3bar ) + PmfHDL[6,NU_TAU_BAR] = ( np.real(HV[NU_MU_BAR,NU_MU_BAR]) - k3bar ) / ( k4bar - k3bar ) + PmfHDL[7,NU_MU_BAR] = ( np.real(HV[NU_TAU_BAR,NU_TAU_BAR]) - k4bar ) / ( k3bar - k4bar ) + PmfHDL[7,NU_TAU_BAR] = ( np.real(HV[NU_MU_BAR,NU_MU_BAR]) - k4bar ) / ( k3bar - k4bar ) + + + """ The IMO case. Matter state 4 is the electron neutrino, matter state 3bar is the electron + antineutrino. Matter state 2 is the sterile neutrino, matter state 1bar is the sterile + antineutrino. + """ + if self.mix_params.mass_order == MassHierarchy.INVERTED: + k1 = T + np.sqrt(D) + k3 = T - np.sqrt(D) + k2bar = Tbar - np.sqrt(Dbar) + k4bar = Tbar + np.sqrt(Dbar) + + PmfHDL[0,NU_MU] = ( np.real(HV[NU_TAU,NU_TAU]) - k1 ) / ( k3 - k1 ) + PmfHDL[0,NU_TAU] = ( np.real(HV[NU_MU,NU_MU]) - k1 ) / ( k3 - k1 ) + PmfHDL[1,NU_S] = 1 + PmfHDL[2,NU_MU] = ( np.real(HV[NU_TAU,NU_TAU]) - k3 ) / ( k1 - k3 ) + PmfHDL[2,NU_TAU] = ( np.real(HV[NU_MU,NU_MU]) - k3 ) / ( k1 - k3 ) + PmfHDL[3,NU_E] = 1 + + PmfHDL[4,NU_S_BAR] = 1 + PmfHDL[5,NU_MU_BAR] = ( np.real(HV[NU_TAU_BAR,NU_TAU_BAR]) - k2bar ) / ( k4bar - k2bar ) + PmfHDL[5,NU_TAU_BAR] = ( np.real(HV[NU_MU_BAR,NU_MU_BAR]) - k2bar ) / ( k4bar - k2bar ) + PmfHDL[6,NU_E_BAR] = 1 + PmfHDL[7,NU_TAU_BAR] = ( np.real(HV[NU_MU_BAR,NU_MU_BAR]) - k4bar ) / ( k2bar - k4bar ) + PmfHDL[7,NU_MU_BAR] = ( np.real(HV[NU_TAU_BAR,NU_TAU_BAR]) - k4bar ) / ( k2bar - k4bar ) + + return PmfHDL + +############################################################################### + +class NoTransformation(FlavorTransformation): + """Survival probabilities for no oscillation case.""" + + def __init__(self): + pass + + def __str__(self): + return f'NoTransformation' + + def get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities. + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + p : 6 x 6 array or array of 6 x 6 arrays + """ + p = np.zeros((6,6,len(E))) + for f in ThreeFlavor: + p[f,f] = np.ones(len(E)) + + return p + +############################################################################### + +class CompleteExchange(FlavorTransformation): + """Survival probabilities for the case when the electron flavors + are half exchanged with the mu flavors and the half with the tau flavors. + """ + + def __init__(self): + pass + + def __str__(self): + return f'CompleteExchange' + + def get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities. + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + p : 6 x 6 array or array of 6 x 6 arrays + """ + p = np.zeros((6,6,len(E))) + + p[ThreeFlavor.NU_E,ThreeFlavor.NU_MU] = np.ones(len(E)) / 2 + p[ThreeFlavor.NU_E,ThreeFlavor.NU_TAU] = np.ones(len(E)) / 2 + p[ThreeFlavor.NU_MU,ThreeFlavor.NU_E] = np.ones(len(E)) / 2 + p[ThreeFlavor.NU_MU,ThreeFlavor.NU_MU] = np.ones(len(E)) / 2 + p[ThreeFlavor.NU_TAU,ThreeFlavor.NU_E] = np.ones(len(E)) / 2 + p[ThreeFlavor.NU_TAU,ThreeFlavor.NU_TAU] = np.ones(len(E)) / 2 + + p[ThreeFlavor.NU_E_BAR,ThreeFlavor.NU_MU_BAR] = np.ones(len(E)) / 2 + p[ThreeFlavor.NU_E_BAR,ThreeFlavor.NU_TAU_BAR] = np.ones(len(E)) / 2 + p[ThreeFlavor.NU_MU_BAR,ThreeFlavor.NU_E_BAR] = np.ones(len(E)) / 2 + p[ThreeFlavor.NU_MU_BAR,ThreeFlavor.NU_MU_BAR] = np.ones(len(E)) / 2 + p[ThreeFlavor.NU_TAU_BAR,ThreeFlavor.NU_E_BAR] = np.ones(len(E)) / 2 + p[ThreeFlavor.NU_TAU_BAR,ThreeFlavor.NU_TAU_BAR] = np.ones(len(E)) / 2 + + return p + +############################################################################### + +class AdiabaticMSW(ThreeFlavorTransformation): + """Adiabatic MSW effect.""" + + def __init__(self, mix_params = None): + """Initialize flavor transformation + + Parameters + ---------- + mix_params : ThreeFlavorMixingParameters instance or None + """ + super().__init__(mix_params) + + def __str__(self): + return f'AdiabaticMSW_' + str(self.mix_params.mass_order) + + def get_SNprobabilities(self, t, E): + """neutrino and antineutrino transition probabilities in the supernova + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + Pmf : array of 6 x 6 arrays + """ + PHDL = self.Pmf_HighDensityLimit() + Pmf = np.empty((6,6,len(E))) + + for m in range(len(E)): + Pmf[:,:,m] = PHDL[:,:] + + return Pmf + + def get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities from the + neutrinosphere to Earth + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + p : array of 6 x 6 arrays + """ + Pmf = self.get_SNprobabilities(t,E) + D = ThreeFlavorNoEarthMatter(self.mix_params).get_probabilities(t,E) + + p = np.empty((6,6,len(E))) + for m in range(len(E)): + p[:,:,m] = D[:,:,m] @ Pmf[:,:,m] + + return p + +############################################################################### + +class NonAdiabaticMSWH(ThreeFlavorTransformation): + """Nonadiabatic MSW H resonance""" + + def __init__(self, mix_params = None): + """Initialize flavor transformation + + Parameters + ---------- + mix_params : ThreeFlavorMixingParameters instance or None + """ + super().__init__(mix_params) + + def __str__(self): + return f'NonAdiabaticMSWH_' + str(self.mix_params.mass_order) + + def get_SNprobabilities(self, t, E): + """neutrino and antineutrino transition probabilities in the supernova + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + Pmf : 6 x 6 array + """ + PHDL = self.Pmf_HighDensityLimit() + Pmf = np.empty((r,6,len(E))) + + for m in range(len(E)): + Pmf[:,:,m] = PHDL[:,:] + if self.mix_params.mass_order == MassHierarchy.NORMAL: + for f in ThreeFlavor: + Pmf[1,f,m], Pmf[2,f,m] = Pmf[2,f,m], Pmf[1,f,m] + if self.mix_params.mass_order == MassHierarchy.INVERTED: + for f in ThreeFlavor: + Pmf[3,f,m], Pmf[5,f,m] = Pmf[5,f,m], Pmf[3,f,m] + + return Pmf + + def get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities. + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + p : 6 x 6 array or array of 6 x 6 arrays + """ + + Pmf = self.get_SNprobabilities(t,E) + D = ThreeFlavorNoEarthMatter(self.mix_params).get_probabilities(t,E) + + p = np.empty((6,6,len(E))) + for m in range(len(E)): + p[:,:,m] = D[:,:,m] @ Pmf[:,:,m] + + return p + +############################################################################### + +class TwoFlavorDecoherence(ThreeFlavorTransformation): + """equal mixing of whatever two matter states form the MSW H resonance""" + + def __init__(self, mix_params = None): + """Initialize flavor transformation + + Parameters + ---------- + mix_params : ThreeFlavorMixingParameters instance or None + """ + super().__init__(mix_params) + + def __str__(self): + return f'TwoFlavorDecoherence_' + str(self.mix_params.mass_order) + + def get_SNprobabilities(self, t, E): + """neutrino and antineutrino transition probabilities. + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + p : 6 x 6 array or array of 6 x 6 arrays + """ + PHDL = self.Pmf_HighDensityLimit(E) + Pmf = np.empty((6,6,len(E))) + + for m in range(len(E)): + Pmf[:,:,m] = PHDL[:,:] + if self.mix_params.mass_order == MassHierarchy.NORMAL: + for f in ThreeFlavor: + Pmf[1,f,m], Pmf[2,f,m] = 0.5 * ( Pmf[1,f,m] + Pmf[2,f,m] ), 0.5 * ( Pmf[1,f,m] + Pmf[2,f,m] ) + if self.mix_params.mass_order == MassHierarchy.INVERTED: + for f in ThreeFlavor: + Pmf[3,f,m], Pmf[5,f,m] = 0.5 * ( Pmf[3,f,m] + Pmf[5,f,m] ), 0.5 * ( Pmf[3,f,m] + Pmf[5,f,m] ) + + return Pmf + + def get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities. + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + p : 6 x 6 array or array of 6 x 6 arrays + """ + + Pmf = self.get_SNprobabilities(t,E) + D = ThreeFlavorNoEarthMatter(self.mix_params).get_probabilities(t,E) + + p = np.empty((6,6,len(E))) + for m in range(len(E)): + p[:,:,m] = D[:,:,m] @ Pmf[:,:,m] + + return p + +############################################################################### + +class ThreeFlavorDecoherence(ThreeFlavorTransformation): + """Equal mixing of all threen eutrino matter states and antineutrino matter states""" + + def __init__(self): + """Initialize ThreeFlavorTransformation to default case""" + super().__init__(None) + + def __str__(self): + return f'ThreeFlavorDecoherence' + + def get_SNprobabilities(self, t, E): + """neutrino and antineutrino transition probabilities. + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + p : 6 x 6 array or array of 6 x 6 arrays + """ + p = np.ones((6,6,len(E))) / 3 + + return p + + def get_probabilities(self, t, E): + """Equal mixing so Earth matter has no effect""" + return self.get_SNprobabilities(t,E) + +############################################################################### + +try: + import SNOSHEWS +except: + SNOSHEWS = None + +class SNprofile: + """A placeholder class for the density and electron fraction profiles. Currently SNOSHEWS + reads the profiles from files but that might change in the future. + """ + + def __init__(self, rhofilename, Yefilename): + self.rhofilename = rhofilename + self.Yefilename = Yefilename + +class MSWEffect(ThreeFlavorTransformation): + """The MSW effect using a density profile and electron + fraction provided by the user. Uses the SNOSHEWS module. + """ + + def __init__(self, SNprofile, mix_params = None, rmin = None, rmax = None): + """Initialize flavor transformation + + Parameters + ---------- + SNprofile : instance of profile class + mix_params : ThreeFlavorMixingParameters instance or None + rmin : starting radius for calculation. rmin will be corrected by SNOSHEWS to the minimum radius of the profile if rmin is less than that value + rmax : the ending radius of the calculation. rmax will be corrected by SNOSHEWS to the maximum radius of the profile if rmax is greater than that value + """ + self.SNprofile = SNprofile + + super().__init__(mix_params) + + if rmin == None: + rmin = 0 + else: + self.rmin = rmin + + if rmax == None: + rmax = 1e99 + else: + self.rmax = rmax + + super().__init__(mix_params) + + def __str__(self): + return f'MSW_' + str(self.mix_params.mass_order) + + def get_SNprobabilities(self, t, E): + """neutrino and antineutrino transition probabilities. + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + Pmf : array of 6 x 6 matrices + """ + + if SNOSHEWS == None: + print("The SNOSHEWS module cannot be found. Returning results using AdiabaticMSW prescription") + return AdiabticMSW(self.mix_params).get_probabilities(t,E) + + #input data object for SNOSHEWS + ID = SNOSHEWS.InputDataSNOSHEWS() + + ID.outputfilenamestem = "./out/SNOSHEWS" # stem of output filenames + + ID.rmin = self.rmin + ID.rmax = self.rmax + + ID.densityprofile = self.SNprofile.rhofilename # mass density profile + ID.electronfraction = self.SNprofile.Yefilename # electron fraction profile + + ID.NE = len(E) # number of energy bins + E = E.to_value('MeV') + ID.Emin = E[0] # in MeV + ID.Emax = E[-1] # in MeV + + #MixingParameters + ID.deltam_21 = self.mix_params.dm21_2.value # in eV^2 + ID.deltam_32 = self.mix_params.dm32_2.value # in eV^2 + ID.theta12 = self.mix_params.theta12.value # in degrees + ID.theta13 = self.mix_params.theta13.value # in degrees + ID.theta23 = self.mix_params.theta23.value # in degrees + ID.deltaCP = self.mix_params.deltaCP.value # in degrees + + ID.accuracy = 1.01E-009 # controls accuracy of integrator: smaller is more accurate + ID.stepcounterlimit = 10000 # output frequency if outputflag = True: larger is less frequent + ID.outputflag = False # set to True if output is desired + + # Do the calculation. The return is a four dimensional + # array of transition probabilities nu_alpha -> nu_i: + # the index order is matter/antimatter, energy, i, alpha + pSN = SNOSHEWS.Run(ID) + + # restructure the results + Pmf = np.zeros((6,6,ID.NE)) + + for m in range(ID.NE): + Pmf[0,ThreeFlavor.NU_E,m] = pSN[0][m][0][0] + Pmf[1,ThreeFlavor.NU_E,m] = pSN[0][m][1][0] + Pmf[2,ThreeFlavor.NU_E,m] = pSN[0][m][2][0] + Pmf[0,ThreeFlavor.NU_MU,m] = pSN[0][m][0][1] + Pmf[1,ThreeFlavor.NU_MU,m] = pSN[0][m][1][1] + Pmf[2,ThreeFlavor.NU_MU,m] = pSN[0][m][2][1] + Pmf[0,ThreeFlavor.NU_TAU,m] = pSN[0][m][0][2] + Pmf[1,ThreeFlavor.NU_TAU,m] = pSN[0][m][1][2] + Pmf[2,ThreeFlavor.NU_TAU,m] = pSN[0][m][2][2] + + Pmf[3,ThreeFlavor.NU_E_BAR,m] = pSN[1][m][0][0] + Pmf[4,ThreeFlavor.NU_E_BAR,m] = pSN[1][m][1][0] + Pmf[5,ThreeFlavor.NU_E_BAR,m] = pSN[1][m][2][0] + Pmf[3,ThreeFlavor.NU_MU_BAR,m] = pSN[1][m][0][1] + Pmf[4,ThreeFlavor.NU_MU_BAR,m] = pSN[1][m][1][1] + Pmf[5,ThreeFlavor.NU_MU_BAR,m] = pSN[1][m][2][1] + Pmf[3,ThreeFlavor.NU_TAU_BAR,m] = pSN[1][m][0][2] + Pmf[4,ThreeFlavor.NU_TAU_BAR,m] = pSN[1][m][1][2] + Pmf[5,ThreeFlavor.NU_TAU_BAR,m] = pSN[1][m][2][2] + + return Pmf + + def get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities. + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + p : 6 x 6 array or array of 6 x 6 arrays + """ + p = np.empty((6,6,len(E))) + + Pmf = self.get_SNprobabilities(t,E) + D = ThreeFlavorNoEarthMatter(self.mix_params).get_probabilities(t,E) + + # multiply the D matrix and the Pmf matrix together + for m in range(len(E)): + p[:,:,m] = self.D[:,:,m] @ Pmf[:,:,m] + + return p + +############################################################################### + +class AdiabaticMSWes(FourFlavorTransformation): + """A four-neutrino mixing prescription. The assumptions used are that: + + 1. the fourth neutrino mass is the heaviest but not so large + that the electron-sterile resonances are inside the neutrinosphere; + 2. the outer electron-sterile MSW resonance is adiabatic; + 3. the inner electron-sterile MSW resonance (where the electron + fraction = 1/3) is non-adiabatic. + + For further insight see, for example, Esmaili, Peres, and Serpico, + Phys. Rev. D 90, 033013 (2014). + """ + def __init__(self, mix_params = None): + """Initialize flavor transformation + + Parameters + ---------- + mix_params : FourFlavorMixingParameters instance or None + """ + super().__init__(mix_params) + + def __str__(self): + return f'AdiabaticMSWes_' + str(self.mix_params.mass_order) + + def get_SNprobabilities(self, t, E): + """neutrino and antineutrino transition probabilities. + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + Pmf : 8 x 8 matrix + """ + PHDL = self.Pmf_HighDensityLimit() + Pmf = np.empty((8,8,len(E))) + + for m in range(len(E)): + Pmf[:,:,m] = PHDL[:,:] + + return Pmf + + + def get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities. + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + p : 6 x 6 matrix + """ + p = np.empty((8,8,len(E))) + + D = FourFlavorNoEarthMatter(self.mix_params).get_probabilities(t,E) + Pmf = self.get_SNprobabilities(t,E) + + # multiply the D matrix and the Pmf matrix together + for m in range(len(E)): + p[:,:,m] = D[:,:,m] @ Pmf[:,:,m] + + # remove sterile rows/columns: A is a 6 x 8 matrix + A = Remove_Steriles() + p = A @ p @ np.transpose(A) + + return p + +############################################################################### + +class NonAdiabaticMSWes(FourFlavorTransformation): + """A four-neutrino mixing prescription. The assumptions used are that: + + 1. the fourth neutrino mass is the heaviest but not so large + that the electron-sterile resonances are inside the neutrinosphere; + 2. the outer electron-sterile MSW resonance is non-adiabatic; + 3. the inner electron-sterile MSW resonance (where the electron + fraction = 1/3) is non-adiabatic. + + For further insight see, for example, Esmaili, Peres, and Serpico, + Phys. Rev. D 90, 033013 (2014). + """ + def __init__(self, mix_params = None): + """Initialize flavor transformation + + Parameters + ---------- + mix_params : FourFlavorMixingParameters instance or None + """ + super().__init__(mix_params) + + def __str__(self): + return f'NonAdiabaticMSWes_' + str(self.mix_params.mass_order) + + def get_SNprobabilities(self, t, E): + """neutrino and antineutrino transition probabilities. + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + Pmf : an 8 x 8 matrix + """ + PHDL = self.Pmf_HighDensityLimit() + Pmf = np.empty((8,8,len(E))) + + for m in range(len(E)): + Pmf[:,:,m] = PHDL[:,:] + + if self.mix_params.mass_order == MassHierarchy.NORMAL: + for f in ThreeFlavor: + Pmf[2,f,m], Pmf[3,f,m] = Pmf[3,f,m], Pmf[2,f,m] + Pmf[5,f,m], Pmf[6,f,m], Pmf[7,f,m] = Pmf[6,f,m], Pmf[7,f,m], Pmf[5,f,m] + + if self.mix_params.mass_order == MassHierarchy.INVERTED: + for f in ThreeFlavor: + Pmf[1,f,m], Pmf[3,f,m] = Pmf[3,f,m], Pmf[1,f,m] + Pmf[4,f,m], Pmf[5,f,m], Pmf[7,f,m] = Pmf[5,f,m], Pmf[7,f,m], Pmf[4,f,m] + + return Pmf + + def get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities. + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + p : 6 x 6 array or array of 6 x 6 arrays + """ + + D = FourFlavorNoEarthMatter(self.mix_params).get_probabilities(t,E) + Pmf = self.get_SNprobabilities(t,E) + + p = np.empty((8,8,len(E))) + for m in range(len(E)): + p[:,:,m] = D[:,:,m] @ Pmf[:,:,m] + + # remove sterile rows/columns + A = Remove_Steriles() + p = A @ p @ np.transpose(A) + + return p + +############################################################################### + +""" The modifier transformtions. These cannot be used by themselves but can be + combined with a flavor transformtion from those defined above using the Catenate class +""" + +class NeutrinoDecay(ThreeFlavorTransformation): + """Decay effect, where the heaviest neutrino decays to the lightest + neutrino. For a description and typical parameters, see A. de Gouvêa et al., + PRD 101:043013, 2020, arXiv:1910.01127. + """ + def __init__(self, mix_params=None, mass=1*u.eV/c.c**2, tau=1*u.day, dist=10*u.kpc): + """Initialize transformation matrix. + + Parameters + ---------- + mix_params : ThreeFlavorMixingParameters instance or None + mass : astropy.units.quantity.Quantity + Mass of the heaviest neutrino; expect in eV/c^2. + tau : astropy.units.quantity.Quantity + Lifetime of the heaviest neutrino. + dist : astropy.units.quantity.Quantity + Distance to the supernova. + """ + super().__init__(mix_params) + + self.m = mass + self.tau = tau + self.d = dist + + def __str__(self): + return f'NeutrinoDecay_' + str(self.mix_params.mass_order) + + def gamma(self, E): + """Decay width of the heaviest neutrino mass. + + Parameters + ---------- + E : float + Energy of the nu3. + + Returns + ------- + Gamma : float + Decay width of the neutrino mass, in units of 1/length. + + :meta private: + """ + return self.m*c.c / (E*self.tau) + + def get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities. + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + PND : an array of 6 x 6 matrices + """ + decay_factor = np.exp(-self.gamma(E)*self.d) + PND = np.zeros((6,6,len(E))) + + if self.mix_params.mass_order == MassHierarchy.NORMAL: + PND[0,0] = 1 + PND[0,2] = 1 - decay_factor + PND[1,1] = 1 + PND[2,2] = decay_factor + + if self.mix_params.mass_order == MassHierarchy.INVERTED: + PND[0,0] = 1 + PND[1,1] = decay_factor + PND[2,1] = 1 - decay_factor + PND[2,2] = 1 + + for i in range(3): + for j in range(3): + PND[i+3,j+3] = PND[i,j] + + return PND + +############################################################################### + +class QuantumDecoherence(ThreeFlavorTransformation): + """Quantum Decoherence, where propagation in vacuum leads to equipartition + of states. For a description and typical parameters, see M. V. dos Santos et al., + 2023, arXiv:2306.17591. + """ + def __init__(self, mix_params=None, Gamma3=1e-27*u.eV, Gamma8=1e-27*u.eV, dist=10*u.kpc, n=0, E0=10*u.MeV): + """Initialize transformation matrix. + + Parameters + ---------- + mix_params : ThreeFlavorMixingParameters instance or None + + Gamma3 : astropy.units.quantity.Quantity + Quantum decoherence parameter; expect in eV. + Gamma8 : astropy.units.quantity.Quantity + Quantum decoherence parameter; expect in eV. + dist : astropy.units.quantity.Quantity + Distance to the supernova. + n : float + Exponent of power law for energy dependent quantum decoherence parameters, + i.e. Gamma = Gamma0*(E/E0)**n. If not specified, it is taken as zero. + E0 : astropy.units.quantity.Quantity + Reference energy in the power law Gamma = Gamma0*(E/E0)**n. If not specified, + it is taken as 10 MeV. Note that if n = 0, quantum decoherence parameters are independent + of E0. + """ + super().__init__(mix_params) + + self.Gamma3 = (Gamma3 / (c.hbar.to('eV s') * c.c)).to('1/kpc') + self.Gamma8 = (Gamma8 / (c.hbar.to('eV s') * c.c)).to('1/kpc') + self.d = dist + self.n = n + self.E0 = E0 + + def __str__(self): + return f'QuantumDecoherence_' + str(self.mix_params.mass_order) + + def get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities. + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + PQD : an array of 6 x 6 matrices + """ + PQD = np.zeros((6,6,len(E))) + + PQD[1,1] = 1/3 + 1/2 * np.exp(-(self.Gamma3 * (E/self.E0)**self.n + self.Gamma8 * (E/self.E0)**self.n / 3) * self.d) \ + + 1/6 * np.exp(-self.Gamma8 * (E/self.E0)**self.n * self.d) + + PQD[1,2] = 1/3 - 1/2 * np.exp(-(self.Gamma3 * (E/self.E0)**self.n + self.Gamma8 * (E/self.E0)**self.n / 3) * self.d) \ + + 1/6 * np.exp(-self.Gamma8 * (E/self.E0)**self.n * self.d) + + PQD[1,3] = 1/3 - 1/3 * np.exp(-self.Gamma8 * (E/self.E0)**self.n * self.d) + + PQD[2,1] = PQD[1,2] + PQD[2,2] = PQD[1,1] + PQD[2,3] = PQD[1,3] + + PQD[3,1] = PQD[1,3] + PQD[3,2] = PQD[2,3] + + PQD[3,3] = 1/3 + 2/3 * np.exp(-self.Gamma8 * (E/self.E0)**self.n * self.d) + + for i in range(3): + for j in range(3): + PQD[i+3,j+3] = PQD[i,j] + + return PQD + +############################################################################### + +class ThreeFlavorNoEarthMatter(ThreeFlavorTransformation): + + def __init__(self, mix_params = None ): + """Initialize flavor transformation + + Parameters + ---------- + mix_params : ThreeFlavorMixingParameters instance or None + """ + super().__init__(mix_params) + + def __str__(self): + pass + + def get_probabilities(self, t, E): + """the D matrix for the case of no Earth matter effects and three neutrino flavors + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + D : an array of length of the E array of 6 x 6 matrices + """ + U = self.mix_params.VacuumMixingMatrix() + + D = np.zeros((6,6,len(E))) # note the first index is a flavor, the second is a mass state + for f in ThreeFlavor: + for i in range(6): + for m in range(len(E)): + D[f,i,m] = float(np.abs(U[f,i])**2) + + return D + +############################################################################### + +class FourFlavorNoEarthMatter(FourFlavorTransformation): + + def __init__(self, mix_params = None ): + """Initialize flavor transformation + + Parameters + ---------- + mix_params : FourFlavorMixingParameters instance or None + """ + super().__init__(mix_params) + + def __str__(self): + pass + + def get_probabilities(self, t, E): + """the D matrix for the case of no Earth matter effects and four neutrino flavors + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + D : an array of 8 x 8 matrices of length equal to the length of the E array + """ + U = self.mix_params.VacuumMixingMatrix() + + D = np.zeros((8,8,len(E))) # note the first index is a flavor, the second is a mass state + for f in FourFlavor: + for i in range(8): + for m in range(len(E)): + D[f,i,m] = float(np.abs(U[f,i])**2) + + return D + +############################################################################### + +try: + import EMEWS +except: + EMEWS = None + +class EarthMatter(ThreeFlavorTransformation): + + def __init__(self, SNAltAz, mix_params = None ): + """Initialize flavor transformation + + Parameters + ---------- + mix_params : ThreeFlavorMixingParameters instance or None + SNAltAz : astropy AltAz object + """ + super().__init__(mix_params) + + self.SNAltAz = SNAltAz + + self.prior_E = None # used to store energy array from previous calls to get_probabilities + self.prior_D = None + + def __str__(self): + return f'EarthMatter_' + str(self.mix_params.mass_order) + + def get_probabilities(self, t, E): + """the D matrix for the case of Earth matter effects + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + D : an array of length of the E array with each element being a 6 x 6 matrix + """ + if EMEWS == None: + print("The EMEWS module cannot be found. Results do not include the Earth-matter effect.") + return ThreeFlavorNoEarthMatter(self.mix_params).get_probabilities(t,E) + + if self.prior_E != None: + if u.allclose(self.prior_E, E) == True: + return self.prior_D + + self.prior_D = np.zeros((6,6,len(E))) + self.prior_E = E + + #input data object for EMEWS + ID = EMEWS.InputDataEMEWS() + + ID.altitude = self.SNAltAz.alt.deg + ID.azimuth = self.SNAltAz.az.deg + + ID.outputfilenamestem = "./out/EMEWS:PREM" # stem of output filenames + ID.densityprofile = "./PREM.rho.dat" # PREM density profile + ID.electronfraction = "./PREM.Ye.dat" # Electron fraction of density profile + + ID.NE = len(E) # number of energy bins + E = E.to_value('MeV') + ID.Emin = E[0] # in MeV + ID.Emax = E[-1] # in MeV + + #MixingParameters + ID.deltam_21 = self.mix_params.dm21_2.value # in eV^2 + ID.deltam_32 = self.mix_params.dm32_2.value # in eV^2 + ID.theta12 = self.mix_params.theta12.value # in degrees + ID.theta13 = self.mix_params.theta13.value # in degrees + ID.theta23 = self.mix_params.theta23.value # in degrees + ID.deltaCP = self.mix_params.deltaCP.value # in degrees + + ID.accuracy = 1.01E-007 # controls accuracy of integrtaor: smaller is more accurate + ID.stepcounterlimit = 1 # output frequency if outputflag = True: larger is less frequent + ID.outputflag = False # set to True if output is desired + + #matrix from EMEWS needs to be rearranged to match SNEWPY flavor indicii ordering + Pfm = EMEWS.Run(ID) + + for m in range(len(E)): + self.prior_D[ThreeFlavor.NU_E,0,m] = Pfm[0][m][0][0] + self.prior_D[ThreeFlavor.NU_E,1,m] = Pfm[0][m][0][1] + self.prior_D[ThreeFlavor.NU_E,2,m] = Pfm[0][m][0][2] + self.prior_D[ThreeFlavor.NU_MU,0,m] = Pfm[0][m][1][0] + self.prior_D[ThreeFlavor.NU_MU,1,m] = Pfm[0][m][1][1] + self.prior_D[ThreeFlavor.NU_MU,2,m] = Pfm[0][m][1][2] + self.prior_D[ThreeFlavor.NU_TAU,0,m] = Pfm[0][m][2][0] + self.prior_D[ThreeFlavor.NU_TAU,1,m] = Pfm[0][m][2][1] + self.prior_D[ThreeFlavor.NU_TAU,2,m] = Pfm[0][m][2][2] + + self.prior_D[ThreeFlavor.NU_E_BAR,3,m] = Pfm[1][m][0][0] + self.prior_D[ThreeFlavor.NU_E_BAR,4,m] = Pfm[1][m][0][1] + self.prior_D[ThreeFlavor.NU_E_BAR,5,m] = Pfm[1][m][0][2] + self.prior_D[ThreeFlavor.NU_MU_BAR,3,m] = Pfm[1][m][1][0] + self.prior_D[ThreeFlavor.NU_MU_BAR,4,m] = Pfm[1][m][1][1] + self.prior_D[ThreeFlavor.NU_MU_BAR,5,m] = Pfm[1][m][1][2] + self.prior_D[ThreeFlavor.NU_TAU_BAR,3,m] = Pfm[1][m][2][0] + self.prior_D[ThreeFlavor.NU_TAU_BAR,4,m] = Pfm[1][m][2][1] + self.prior_D[ThreeFlavor.NU_TAU_BAR,5,m] = Pfm[1][m][2][2] + + return self.prior_D + +############################################################################### + +class Catenate: + """Catenate flavor transformation effects together.""" + + def __init__(self, SNTransformation, InVacuumTransformation = None, AtEarthTransformation = None): + """ + Parameters + ---------- + SNTransformation : a FlavorTransformation object that describes the transformation that occurs in the supernova + InVacuumTransformation : a FlavorTransformation object that describes the transforamtion that occurs in the vacuum + AtEarthTransformation : a FlavorTransformation object that describes the transformation that occurs at Earth + The order is that SNTransformation is applied first, then InVacuumTransformation, then AtEarthTransformation + """ + self.transform1 = SNTransformation + + if InVacuumTransformation == None: + self.transform2 = NoTransformation() + else: + self.transform2 = InVacuumTransformation + + if AtEarthTransformation == None: + self.transform3 = ThreeFlavorNoEarthMatter(self.transform1.mix_params) + else: + self.transform3 = AtEarthTransformation + + + def __str__(self): + return str(self.transform1) + '+' + str(self.transform2) + '+' + str(self.transform3) + + def get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities. + + Parameters + ---------- + t : float or ndarray + List of times. + E : float or ndarray + List of energies. + + Returns + ------- + p : 6 x 6 array or array of 6 x 6 arrays + """ + p1 = self.transform1.get_SNprobabilities(t,E) + p2 = self.transform2.get_probabilities(t,E) + p3 = self.transform3.get_probabilities(t,E) + + p = np.empty((6,6,len(E))) + for m in range(len(E)): + p[:,:,m] = p3[:,:,m] @ p2[:,:,m] @ p1[:,:,m] + + return p From bd923c7c80b4fe06b2541a566414389dd4afe5c7 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Fri, 7 Jun 2024 14:56:00 +0300 Subject: [PATCH 02/42] Add missing imports, remove Two/Three/FourFlavor definitions from neutrino.py --- python/snewpy/flavor_transformation.py | 23 +--- python/snewpy/neutrino.py | 145 +++++++++++++++++++++---- 2 files changed, 132 insertions(+), 36 deletions(-) diff --git a/python/snewpy/flavor_transformation.py b/python/snewpy/flavor_transformation.py index e5e6802ab..1fba0dbc9 100644 --- a/python/snewpy/flavor_transformation.py +++ b/python/snewpy/flavor_transformation.py @@ -12,7 +12,8 @@ from astropy import constants as c from astropy.coordinates import AltAz -from .neutrino import MassHierarchy, ThreeFlavor, FourFlavor +from .neutrino import MassHierarchy +from .flavor import ThreeFlavor, FourFlavor, FlavorMatrix ############################################################################### @@ -20,7 +21,7 @@ class FlavorTransformation(ABC): """Generic interface to compute neutrino and antineutrino survival probability.""" @abstractmethod - def get_probabilities(self, t, E): + def get_probabilities(self, t, E)->FlavorMatrix: """neutrino and antineutrino transition probabilities. Parameters @@ -32,8 +33,8 @@ def get_probabilities(self, t, E): Returns ------- - p : a N x N array, or an array of N x N arrays - where N is either 6 or 8 + p : a [N x N] array or [N x N x len(E) x len(t)] array + where N is number of neutrino flavors(6 or 8) """ pass @@ -248,19 +249,7 @@ def __str__(self): return f'NoTransformation' def get_probabilities(self, t, E): - """neutrino and antineutrino transition probabilities. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - p : 6 x 6 array or array of 6 x 6 arrays - """ + p = np.zeros((6,6,len(E))) for f in ThreeFlavor: p[f,f] = np.ones(len(E)) diff --git a/python/snewpy/neutrino.py b/python/snewpy/neutrino.py index e1612360f..919846418 100644 --- a/python/snewpy/neutrino.py +++ b/python/snewpy/neutrino.py @@ -7,8 +7,7 @@ from typing import Optional import numpy as np from collections.abc import Mapping -from .flavor import TwoFlavor as Flavor - +from .flavor import ThreeFlavor, FourFlavor class MassHierarchy(IntEnum): """Neutrino mass ordering: ``NORMAL`` or ``INVERTED``.""" NORMAL = 1 @@ -23,10 +22,16 @@ def derive_from_dm2(cls, dm12_2, dm32_2, dm31_2): return MassHierarchy.NORMAL else: return MassHierarchy.INVERTED - + + def __str__(self): + if self.value == MassHierarchy.NORMAL: + return f'NMO' + if self.value == MassHierarchy.INVERTED: + return f'IMO' + @dataclass -class MixingParameters3Flavor(Mapping): +class ThreeFlavorMixingParameters(Mapping): """Mixing angles and mass differences, assuming three neutrino flavors. This class contains the default values used throughout SNEWPY, which are based on `NuFIT 5.0 `_ results from July 2020, @@ -110,18 +115,61 @@ def get_mass_squared_differences(self): """ return (self.dm21_2, self.dm31_2, self.dm32_2) + + def VacuumMixingMatrix(self): + """The vacuum mixing matrix given the mixing paramters + N.B. This is a 6 x 6 matrix + """ + + U = self.ComplexRotationMatrix(1,2,self.theta23,0) \ + @ self.ComplexRotationMatrix(0,2,self.theta13,self.deltaCP) \ + @ self.ComplexRotationMatrix(0,1,self.theta12,0) + + """Reorder rows to match SNEWPY's flavor ordering convention""" + U[ThreeFlavor.NU_E], U[ThreeFlavor.NU_MU], U[ThreeFlavor.NU_TAU], \ + U[ThreeFlavor.NU_E_BAR], U[ThreeFlavor.NU_MU_BAR], U[ThreeFlavor.NU_TAU_BAR] = \ + U[0], U[1], U[2], U[3], U[4], U[5] + + return U + + + def ComplexRotationMatrix(self,i,j,theta,phase): + """A complex rotation matrix. N.B. the minus sign in the complex exponential matches PDG convention""" + V = np.zeros((6,6),dtype = 'complex_') + for k in range(6): + V[k,k] = 1 + + V[i,i] = np.cos(theta) + V[j,j] = V[i,i] + V[i,j] = np.sin(theta) * ( np.cos(phase) - 1j * np.sin(phase) ) + V[j,i] = - np.conjugate(V[i,j]) + + V[i+3,i+3] = V[i,i] + V[j+3,j+3] = V[j,j] + V[i+3,j+3] = np.conjugate(V[i,j]) + V[j+3,i+3] = np.conjugate(V[j,i]) + + return V + + @dataclass -class MixingParameters4Flavor(MixingParameters3Flavor): +class FourFlavorMixingParameters(ThreeFlavorMixingParameters): """A class for four flavor neutrino mixing. - ..Note: it is an extension of :class:`MixingParameters3Flavor`, and can be constructed using it: + ..Note: it is an extension of :class:`ThreeFlavorMixingParameters`, and can be constructed using it: - >>> pars_3f = MixingParameters() #standard 3flavor mixing - >>> pars_4f = MixingParameters4Flavor(**pars_3f, theta14=90<>> pars_3f = ThreeFlavorMixingParameters() #standard 3flavor mixing + >>> pars_4f = FpourFlavorMixingParameters(**pars_3f, theta14=90< Date: Fri, 7 Jun 2024 16:39:51 +0300 Subject: [PATCH 03/42] Adding 'zeros' and 'T' methods for FlavorMatrix --- python/snewpy/flavor.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/python/snewpy/flavor.py b/python/snewpy/flavor.py index 9ed92eaaf..f371a6120 100644 --- a/python/snewpy/flavor.py +++ b/python/snewpy/flavor.py @@ -131,6 +131,10 @@ def shape(self): @property def flavor(self): return self.flavor_out + @property + def T(self): + "transposed version of the matrix: reverse the flavor dimensions" + return FlavorMatrix(array = self.swapaxes(0,1), flavor = self.flavor_in, from_flavor=self.flavor_out) @classmethod def eye(cls, flavor:FlavorScheme, from_flavor:FlavorScheme = None): @@ -138,7 +142,12 @@ def eye(cls, flavor:FlavorScheme, from_flavor:FlavorScheme = None): shape = (len(from_flavor), len(flavor)) data = np.eye(*shape) return cls(data, flavor, from_flavor) - + @classmethod + def zeros(cls, flavor:FlavorScheme, from_flavor:FlavorScheme = None): + from_flavor = from_flavor or flavor + shape = (len(from_flavor), len(flavor)) + data = np.zeros(shape) + return cls(data, flavor, from_flavor) @classmethod def from_function(cls, flavor:FlavorScheme, from_flavor:FlavorScheme = None): """A decorator for creating the flavor matrix from the given function""" From de21f16d384a9a534ef9481ae329cb58e1aaa6db Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Fri, 7 Jun 2024 16:40:40 +0300 Subject: [PATCH 04/42] Removing repeating code, using FlavorMatrices --- python/snewpy/flavor_transformation.py | 163 ++++--------------------- 1 file changed, 26 insertions(+), 137 deletions(-) diff --git a/python/snewpy/flavor_transformation.py b/python/snewpy/flavor_transformation.py index 1fba0dbc9..ad5f52947 100644 --- a/python/snewpy/flavor_transformation.py +++ b/python/snewpy/flavor_transformation.py @@ -13,13 +13,16 @@ from astropy.coordinates import AltAz from .neutrino import MassHierarchy -from .flavor import ThreeFlavor, FourFlavor, FlavorMatrix +from .flavor import TwoFlavor,ThreeFlavor, FourFlavor, FlavorMatrix ############################################################################### class FlavorTransformation(ABC): """Generic interface to compute neutrino and antineutrino survival probability.""" + def __str__(self): + return self.__class__.__name__ + @abstractmethod def get_probabilities(self, t, E)->FlavorMatrix: """neutrino and antineutrino transition probabilities. @@ -55,7 +58,9 @@ def __init__(self, mix_params): else: self.mix_params = mix_params - + def __str__(self): + return self.__class__.__name__+"_"+str(self.mix_params.mass_order) + def Pmf_HighDensityLimit(self): """ The probability that a given flavor state is a particular matter state in the infinite density limit. @@ -148,6 +153,9 @@ def __init__(self, mix_params): else: self.mix_params = mix_params + def __str__(self): + return self.__class__.__name__+"_"+str(self.mix_params.mass_order) + def Pmf_HighDensityLimit(self): """ The probability that a given flavor state is 'detected' as a particular matter state in the infinite density limit. This method assumes the 4th matter state is the @@ -242,18 +250,8 @@ def Pmf_HighDensityLimit(self): class NoTransformation(FlavorTransformation): """Survival probabilities for no oscillation case.""" - def __init__(self): - pass - - def __str__(self): - return f'NoTransformation' - - def get_probabilities(self, t, E): - - p = np.zeros((6,6,len(E))) - for f in ThreeFlavor: - p[f,f] = np.ones(len(E)) - + def get_probabilities(self, t, E): + p = FlavorMatrix.eye(ThreeFlavor) return p ############################################################################### @@ -262,63 +260,20 @@ class CompleteExchange(FlavorTransformation): """Survival probabilities for the case when the electron flavors are half exchanged with the mu flavors and the half with the tau flavors. """ - - def __init__(self): - pass - - def __str__(self): - return f'CompleteExchange' - def get_probabilities(self, t, E): - """neutrino and antineutrino transition probabilities. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - p : 6 x 6 array or array of 6 x 6 arrays - """ - p = np.zeros((6,6,len(E))) - - p[ThreeFlavor.NU_E,ThreeFlavor.NU_MU] = np.ones(len(E)) / 2 - p[ThreeFlavor.NU_E,ThreeFlavor.NU_TAU] = np.ones(len(E)) / 2 - p[ThreeFlavor.NU_MU,ThreeFlavor.NU_E] = np.ones(len(E)) / 2 - p[ThreeFlavor.NU_MU,ThreeFlavor.NU_MU] = np.ones(len(E)) / 2 - p[ThreeFlavor.NU_TAU,ThreeFlavor.NU_E] = np.ones(len(E)) / 2 - p[ThreeFlavor.NU_TAU,ThreeFlavor.NU_TAU] = np.ones(len(E)) / 2 - - p[ThreeFlavor.NU_E_BAR,ThreeFlavor.NU_MU_BAR] = np.ones(len(E)) / 2 - p[ThreeFlavor.NU_E_BAR,ThreeFlavor.NU_TAU_BAR] = np.ones(len(E)) / 2 - p[ThreeFlavor.NU_MU_BAR,ThreeFlavor.NU_E_BAR] = np.ones(len(E)) / 2 - p[ThreeFlavor.NU_MU_BAR,ThreeFlavor.NU_MU_BAR] = np.ones(len(E)) / 2 - p[ThreeFlavor.NU_TAU_BAR,ThreeFlavor.NU_E_BAR] = np.ones(len(E)) / 2 - p[ThreeFlavor.NU_TAU_BAR,ThreeFlavor.NU_TAU_BAR] = np.ones(len(E)) / 2 - - return p + p = FlavorMatrix.zeros(TwoFlavor) + p['NU_X','NU_E'] = 1 + p['NU_E','NU_X'] = 1 + p['NU_E_BAR','NU_X_BAR'] = 1 + p['NU_X_BAR','NU_E_BAR'] = 1 + return (ThreeFlavor< Date: Fri, 7 Jun 2024 16:41:13 +0300 Subject: [PATCH 05/42] Adding default Flavor import --- python/snewpy/neutrino.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/python/snewpy/neutrino.py b/python/snewpy/neutrino.py index 919846418..6c77605ba 100644 --- a/python/snewpy/neutrino.py +++ b/python/snewpy/neutrino.py @@ -8,6 +8,8 @@ import numpy as np from collections.abc import Mapping from .flavor import ThreeFlavor, FourFlavor +from .flavor import TwoFlavor as Flavor + class MassHierarchy(IntEnum): """Neutrino mass ordering: ``NORMAL`` or ``INVERTED``.""" NORMAL = 1 From c1f6eabd9cee89360f343172c6d8795883d06f34 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Fri, 7 Jun 2024 17:27:29 +0300 Subject: [PATCH 06/42] Use matrix multiplication in model --- python/snewpy/models/base.py | 42 ++++-------------------------------- 1 file changed, 4 insertions(+), 38 deletions(-) diff --git a/python/snewpy/models/base.py b/python/snewpy/models/base.py index 0d81efc0c..188906d59 100644 --- a/python/snewpy/models/base.py +++ b/python/snewpy/models/base.py @@ -12,10 +12,11 @@ from snewpy._model_downloader import LocalFileLoader from snewpy.neutrino import Flavor +from snewpy.flavor import TwoFlavor, ThreeFlavor, FourFlavor from snewpy.flavor_transformation import NoTransformation from functools import wraps -from snewpy.flux import Flux +from snewpy.flux import Flux, Container from pathlib import Path def _wrap_init(init, check): @@ -119,14 +120,6 @@ def get_initial_spectra(self, t, E, flavors=Flavor): """ pass - def get_initialspectra(self, *args): - """DO NOT USE! Only for backward compatibility! - - :meta private: - """ - warn("Please use `get_initial_spectra()` instead of `get_initialspectra()`!", FutureWarning) - return self.get_initial_spectra(*args) - def get_transformed_spectra(self, t, E, flavor_xform): """Get neutrino spectra after applying oscillation. @@ -144,25 +137,8 @@ def get_transformed_spectra(self, t, E, flavor_xform): dict Dictionary of transformed spectra, keyed by neutrino flavor. """ - initialspectra = self.get_initial_spectra(t, E) - transformed_spectra = {} - - transformed_spectra[Flavor.NU_E] = \ - flavor_xform.prob_ee(t, E) * initialspectra[Flavor.NU_E] + \ - flavor_xform.prob_ex(t, E) * initialspectra[Flavor.NU_X] - - transformed_spectra[Flavor.NU_X] = \ - flavor_xform.prob_xe(t, E) * initialspectra[Flavor.NU_E] + \ - flavor_xform.prob_xx(t, E) * initialspectra[Flavor.NU_X] - - transformed_spectra[Flavor.NU_E_BAR] = \ - flavor_xform.prob_eebar(t, E) * initialspectra[Flavor.NU_E_BAR] + \ - flavor_xform.prob_exbar(t, E) * initialspectra[Flavor.NU_X_BAR] - - transformed_spectra[Flavor.NU_X_BAR] = \ - flavor_xform.prob_xebar(t, E) * initialspectra[Flavor.NU_E_BAR] + \ - flavor_xform.prob_xxbar(t, E) * initialspectra[Flavor.NU_X_BAR] - + initial_spectra = self.get_initial_spectra(t, E) + transformed_spectra = flavor_xform.get_probabilities(t,E) @ (ThreeFlavor< Date: Fri, 7 Jun 2024 17:28:04 +0300 Subject: [PATCH 07/42] Add matrix multiplication to dict --- python/snewpy/flavor.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/python/snewpy/flavor.py b/python/snewpy/flavor.py index f371a6120..7f8c42fe3 100644 --- a/python/snewpy/flavor.py +++ b/python/snewpy/flavor.py @@ -122,7 +122,13 @@ def __matmul__(self, other): except Exception as e: raise ValueError(f"Cannot multiply {self._repr_short()} by {other._repr_short()}") from e elif hasattr(other, '__rmatmul__'): - return other.__rmatmul__(self) + return other.__rmatmul__(self) + elif isinstance(other, dict): + #try to multiply to the dict[Flavor:flux] + #we assume that it has the same Flavor scheme! + array = np.stack([other[flv] for flv in self.flavor_in]) + result = np.tensordot(self.array, array, axes=[1,0]) + return {flv:result[n] for n,flv in enumerate(self.flavor_out)} raise TypeError(f"Cannot multiply object of {self.__class__} by {other.__class__}") #properties @property From ad7674395af4f135993c990abc9b53158a4e3bd0 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Sun, 16 Jun 2024 21:40:14 +0300 Subject: [PATCH 08/42] Using einsum to do the proper Matrix@Flux multiplication --- python/snewpy/flux.py | 12 ++++++++++-- python/snewpy/utils.py | 7 +++++++ 2 files changed, 17 insertions(+), 2 deletions(-) create mode 100644 python/snewpy/utils.py diff --git a/python/snewpy/flux.py b/python/snewpy/flux.py index fc9469153..d62f8ae31 100644 --- a/python/snewpy/flux.py +++ b/python/snewpy/flux.py @@ -62,7 +62,7 @@ from scipy.interpolate import interp1d from enum import IntEnum from functools import wraps - +import snewpy.utils #list of units which will be used as units for decomposition inside the Container snewpy_unit_bases = [u.MeV, u.m, u.s, u.kg] @@ -404,11 +404,19 @@ def __rshift__(self, flavor:FlavorScheme): return self.convert_to_flavor(flavor) def __rmatmul__(self, matrix:FlavorMatrix): + """Multiply this flux by a FlavorMatrix""" if not self._is_full_flavor(): raise RuntimeError(f"Cannot multiply flavor matrix object {self}, expected {len(self.flavor_scheme)} flavors") if matrix.flavor_in!=self.flavor_scheme: raise ValueError(f"Cannot multiply flavor matrix {matrix} by {self} - flavor scheme mismatch!") - array = np.tensordot(matrix.array,self.array, axes=[1,0]) + #apply the multiplication: + #first add the missing dimensions for the matrix (if needed) + f = self.array + m = snewpy.utils.expand_dimensions_to(matrix.array, + ndim=f.ndim+1) + + #do the multiplication + array = np.einsum('ij...,j...->i...',m,f) return Container(array, flavor=matrix.flavor_out, time=self.time, energy=self.energy) class Container(_ContainerBase): #a dictionary holding classes for each unit diff --git a/python/snewpy/utils.py b/python/snewpy/utils.py new file mode 100644 index 000000000..5d3b488bb --- /dev/null +++ b/python/snewpy/utils.py @@ -0,0 +1,7 @@ +import numpy as np + +def expand_dimensions_to(a:np.ndarray, ndim:int)->np.ndarray: + """Expand the dimensions of the array, adding dimensions of len=1 to the right, + so total dimensions equal to `ndim`""" + new_shape = (list(a.shape)+[1]*ndim)[:ndim] + return a.reshape(new_shape) \ No newline at end of file From b9f60980312c245f8f49b6549a54356d57ee099e Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Sun, 16 Jun 2024 22:13:39 +0300 Subject: [PATCH 09/42] Allow creating FlavorMatrices with extra dimensions --- python/snewpy/flavor.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/python/snewpy/flavor.py b/python/snewpy/flavor.py index 7f8c42fe3..03972a6b8 100644 --- a/python/snewpy/flavor.py +++ b/python/snewpy/flavor.py @@ -90,7 +90,7 @@ def __init__(self, self.flavor_out = flavor self.flavor_in = from_flavor or flavor expected_shape = (len(self.flavor_out), len(self.flavor_in)) - if(self.array.shape != expected_shape): + if(self.array.shape[:2] != expected_shape): raise ValueError(f"FlavorMatrix array shape {self.array.shape} mismatch expected {expected_shape}") def _convert_index(self, index): @@ -109,7 +109,9 @@ def _repr_short(self): return f'{self.__class__.__name__}:<{self.flavor_in.__name__}->{self.flavor_out.__name__}> shape={self.shape}' def __repr__(self): - s=self._repr_short()+'\n'+repr(self.array) + s=self._repr_short() + if(len(self.shape)==2): + s+='\n'+repr(self.array) return s def __eq__(self,other): return self.flavor_in==other.flavor_in and self.flavor_out==other.flavor_out and np.allclose(self.array,other.array) From 8e257c25c4a266a0619249717b9ecc923dea06fe Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Sun, 16 Jun 2024 22:21:37 +0300 Subject: [PATCH 10/42] Using einsum instead of tensordot --- python/snewpy/flavor.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/python/snewpy/flavor.py b/python/snewpy/flavor.py index 03972a6b8..644fd21b8 100644 --- a/python/snewpy/flavor.py +++ b/python/snewpy/flavor.py @@ -119,8 +119,12 @@ def __eq__(self,other): def __matmul__(self, other): if isinstance(other, FlavorMatrix): try: - data = np.tensordot(self.array, other.array, axes=[1,0]) - return FlavorMatrix(data, self.flavor_out, from_flavor = other.flavor_in) + m0, m1 = self.array, other.array + ndims = max(m0.ndim, m1.ndim) + m0,m1 = [snewpy.utils.expand_dimensions_to(m, ndim=ndims) for m in [m0,m1]] + array = np.einsum('ij...,jk...->ik...',m,f) + np.tensordot(self.array, other.array, axes=[1,0]) + return FlavorMatrix(array, self.flavor_out, from_flavor = other.flavor_in) except Exception as e: raise ValueError(f"Cannot multiply {self._repr_short()} by {other._repr_short()}") from e elif hasattr(other, '__rmatmul__'): From 5c926eda5a851ad527651ddfb118f67c1d48114b Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Sun, 16 Jun 2024 22:35:58 +0300 Subject: [PATCH 11/42] Fix typo --- python/snewpy/flavor.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/python/snewpy/flavor.py b/python/snewpy/flavor.py index 644fd21b8..d2bff3157 100644 --- a/python/snewpy/flavor.py +++ b/python/snewpy/flavor.py @@ -1,6 +1,7 @@ import enum import numpy as np import typing +import snewpy.utils class EnumMeta(enum.EnumMeta): def __getitem__(cls, key): @@ -122,7 +123,7 @@ def __matmul__(self, other): m0, m1 = self.array, other.array ndims = max(m0.ndim, m1.ndim) m0,m1 = [snewpy.utils.expand_dimensions_to(m, ndim=ndims) for m in [m0,m1]] - array = np.einsum('ij...,jk...->ik...',m,f) + array = np.einsum('ij...,jk...->ik...',m0,m1) np.tensordot(self.array, other.array, axes=[1,0]) return FlavorMatrix(array, self.flavor_out, from_flavor = other.flavor_in) except Exception as e: From b1363a205be697aa56a287fc9dce787e53bcdf28 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Mon, 17 Jun 2024 20:16:37 +0300 Subject: [PATCH 12/42] Implement some transformations with the FlavorMatrix --- python/snewpy/flavor_transformation.py | 76 +++++++------------------- 1 file changed, 19 insertions(+), 57 deletions(-) diff --git a/python/snewpy/flavor_transformation.py b/python/snewpy/flavor_transformation.py index ad5f52947..91951dcce 100644 --- a/python/snewpy/flavor_transformation.py +++ b/python/snewpy/flavor_transformation.py @@ -14,7 +14,7 @@ from .neutrino import MassHierarchy from .flavor import TwoFlavor,ThreeFlavor, FourFlavor, FlavorMatrix - +from .neutrino import MixingParameters, ThreeFlavorMixingParameters, FourFlavorMixingParameters ############################################################################### class FlavorTransformation(ABC): @@ -38,9 +38,12 @@ def get_probabilities(self, t, E)->FlavorMatrix: ------- p : a [N x N] array or [N x N x len(E) x len(t)] array where N is number of neutrino flavors(6 or 8) - """ + """ pass - + + def apply(self, flux): + M = self.get_probabilities(flux.time, flux.energy) + return M@flux ############################################################################### class ThreeFlavorTransformation(FlavorTransformation): @@ -53,11 +56,8 @@ def __init__(self, mix_params): ---------- mix_params : ThreeFlavorMixingParameters instance or None """ - if mix_params == None: - self.mix_params = ThreeFlavorMixingParameters(MassHierarchy.NORMAL) - else: - self.mix_params = mix_params - + self.mix_params = mix_params or MixingParameters(MassHierarchy.NORMAL) + def __str__(self): return self.__class__.__name__+"_"+str(self.mix_params.mass_order) @@ -250,10 +250,13 @@ def Pmf_HighDensityLimit(self): class NoTransformation(FlavorTransformation): """Survival probabilities for no oscillation case.""" - def get_probabilities(self, t, E): + def get_probabilities(self, t, E): p = FlavorMatrix.eye(ThreeFlavor) return p + def apply(self, flux): + return flux + ############################################################################### class CompleteExchange(FlavorTransformation): @@ -273,37 +276,9 @@ def get_probabilities(self, t, E): class AdiabaticMSW(ThreeFlavorTransformation): """Adiabatic MSW effect.""" - def _get_SNprobabilities(self, t, E): - """neutrino and antineutrino transition probabilities in the supernova - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - Pmf : array of 6 x 6 arrays - """ + def get_probabilities(self, t, E): PHDL = self.Pmf_HighDensityLimit() - Pmf = np.empty((6,6,len(E))) - - for m in range(len(E)): - Pmf[:,:,m] = PHDL[:,:] - - return Pmf - - def get_probabilities(self, t, E): - Pmf = self._get_SNprobabilities(t,E) - D = ThreeFlavorNoEarthMatter(self.mix_params).get_probabilities(t,E) - - p = np.empty((6,6,len(E))) - for m in range(len(E)): - p[:,:,m] = D[:,:,m] @ Pmf[:,:,m] - - return p + return FlavorMatrix(PHDL,ThreeFlavor) ############################################################################### @@ -382,7 +357,7 @@ def get_SNprobabilities(self, t, E): return Pmf def get_probabilities(self, t, E): - Pmf = self.get_SNprobabilities(t,E) + Pmf = self.get_SNprobabilities(t,E) D = ThreeFlavorNoEarthMatter(self.mix_params).get_probabilities(t,E) p = np.empty((6,6,len(E))) @@ -936,14 +911,8 @@ def get_probabilities(self, t, E): D : an array of length of the E array of 6 x 6 matrices """ U = self.mix_params.VacuumMixingMatrix() - - D = np.zeros((6,6,len(E))) # note the first index is a flavor, the second is a mass state - for f in ThreeFlavor: - for i in range(6): - for m in range(len(E)): - D[f,i,m] = float(np.abs(U[f,i])**2) - - return D + P = np.abs(U**2) + return FlavorMatrix(P, ThreeFlavor) ############################################################################### @@ -976,15 +945,8 @@ def get_probabilities(self, t, E): D : an array of 8 x 8 matrices of length equal to the length of the E array """ U = self.mix_params.VacuumMixingMatrix() - - D = np.zeros((8,8,len(E))) # note the first index is a flavor, the second is a mass state - for f in FourFlavor: - for i in range(8): - for m in range(len(E)): - D[f,i,m] = float(np.abs(U[f,i])**2) - - return D - + P = np.abs(U**2) + return FlavorMatrix(P, FourFlavor) ############################################################################### try: From 4e2be8e6aaa2b0b5c528b4417b489f944308453f Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Mon, 17 Jun 2024 20:55:17 +0300 Subject: [PATCH 13/42] Closes #342: adding the function --- python/snewpy/flux.py | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/python/snewpy/flux.py b/python/snewpy/flux.py index d62f8ae31..9085cc497 100644 --- a/python/snewpy/flux.py +++ b/python/snewpy/flux.py @@ -66,6 +66,8 @@ #list of units which will be used as units for decomposition inside the Container snewpy_unit_bases = [u.MeV, u.m, u.s, u.kg] +import matplotlib.pyplot as plt + class Axes(IntEnum): """Enum to keep the number number of the array dimension for each axis""" flavor=0, #Flavor dimension @@ -442,7 +444,27 @@ def __class_getitem__(cls, args): name = name or f'{cls.__name__}[{unit}]' cls._unit_classes[unit] = type(name,(cls,),{'unit':unit}) return cls._unit_classes[unit] - + + def plot(flux, projection='energy', styles=None, **kwargs): + if projection=='energy': + fP = flux.integrate_or_sum('time') + else: + fP = flux.integrate_or_sum('energy') + x = fP.__dict__[projection] + if isinstance(styles, dict): + styles = styles.get + elif styles==None: + styles = lambda flv: {'ls':'-' if flv.is_neutrino else ':', + 'color':f'C{flv//2:d}'} + for flv in fP.flavor: + style = styles(flv) + style.update(kwargs) + y = fP[flv].array.squeeze() + plt.plot(x,y,label=flv.to_tex(), **style) + plt.legend(ncols=2) + plt.xlabel(f'{projection}, {x.unit._repr_latex_()}') + plt.ylabel(f'{fP.__class__.__name__}, {x.unit._repr_latex_()}') + return #some standard container classes that can be used for Flux = Container['1/(MeV*s*m**2)', "d2FdEdT"] From 2eebeb332875112fd9264cb1e78538754d246f43 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Thu, 25 Jul 2024 12:06:52 +0300 Subject: [PATCH 14/42] Fix the separator strings --- python/snewpy/models/presn_loaders.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/snewpy/models/presn_loaders.py b/python/snewpy/models/presn_loaders.py index c116207bc..721a8742c 100644 --- a/python/snewpy/models/presn_loaders.py +++ b/python/snewpy/models/presn_loaders.py @@ -41,7 +41,7 @@ class Odrzywolek_2010(SupernovaModel): def __init__(self, filename:str): df = pd.read_csv( self.request_file(filename), - sep='\s+', + sep=r'\s+', skiprows=1, usecols=[1,6,7,8], names=["time","a","alpha","b"], @@ -83,7 +83,7 @@ def __init__(self, filename:str): df = pd.read_csv( self.request_file(filename), comment="#", - sep='\s+', + sep=r'\s+', names=["time","Enu",Flavor.NU_E,Flavor.NU_E_BAR,Flavor.NU_MU,Flavor.NU_MU_BAR], usecols=range(6), ) From 21f21f1925df4793fcc4fe6dd318cd5a625cfdc1 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Thu, 25 Jul 2024 12:07:41 +0300 Subject: [PATCH 15/42] Getting rid of 'Optional' type annotations --- python/snewpy/neutrino.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/python/snewpy/neutrino.py b/python/snewpy/neutrino.py index 193d6cba9..cc1b6423d 100644 --- a/python/snewpy/neutrino.py +++ b/python/snewpy/neutrino.py @@ -164,13 +164,13 @@ class FourFlavorMixingParameters(ThreeFlavorMixingParameters): """ #sterile neutrino mixing angles. theta14: u.Quantity[u.deg] = 0< Date: Thu, 25 Jul 2024 16:07:13 +0300 Subject: [PATCH 16/42] Adding check that flavor schemes are compatible in __matmul__ --- python/snewpy/flavor.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/python/snewpy/flavor.py b/python/snewpy/flavor.py index d2bff3157..2f56c588a 100644 --- a/python/snewpy/flavor.py +++ b/python/snewpy/flavor.py @@ -119,6 +119,7 @@ def __eq__(self,other): def __matmul__(self, other): if isinstance(other, FlavorMatrix): + assert self.flavor_in==other.flavor_out, f"Incompatible spaces {self.flavor_in}!={other.flavor_out}" try: m0, m1 = self.array, other.array ndims = max(m0.ndim, m1.ndim) @@ -146,8 +147,16 @@ def flavor(self): return self.flavor_out @property def T(self): + return self.transpose() + + def transpose(self): "transposed version of the matrix: reverse the flavor dimensions" - return FlavorMatrix(array = self.swapaxes(0,1), flavor = self.flavor_in, from_flavor=self.flavor_out) + return FlavorMatrix(array = self.array.swapaxes(0,1), + flavor = self.flavor_in, from_flavor=self.flavor_out) + def conjugate(self): + "apply complex conjugate" + return FlavorMatrix(array = self.array.conjugate(), + flavor = self.flavor_out, from_flavor=self.flavor_in) @classmethod def eye(cls, flavor:FlavorScheme, from_flavor:FlavorScheme = None): From 48eebe0ad63994503bc895206d94f27d2aa28515 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Thu, 25 Jul 2024 17:05:43 +0300 Subject: [PATCH 17/42] Moving Pmf_HDL function to MixingParameters class --- python/snewpy/flavor_transformation.py | 83 ++------------------ python/snewpy/neutrino.py | 102 +++++++++++++++++++------ 2 files changed, 86 insertions(+), 99 deletions(-) diff --git a/python/snewpy/flavor_transformation.py b/python/snewpy/flavor_transformation.py index 91951dcce..d3a62d0ea 100644 --- a/python/snewpy/flavor_transformation.py +++ b/python/snewpy/flavor_transformation.py @@ -22,7 +22,7 @@ class FlavorTransformation(ABC): def __str__(self): return self.__class__.__name__ - + @abstractmethod def get_probabilities(self, t, E)->FlavorMatrix: """neutrino and antineutrino transition probabilities. @@ -44,6 +44,9 @@ def get_probabilities(self, t, E)->FlavorMatrix: def apply(self, flux): M = self.get_probabilities(flux.time, flux.energy) return M@flux + + def P(self, t, E): + return self.get_probabilities(t,E) ############################################################################### class ThreeFlavorTransformation(FlavorTransformation): @@ -60,81 +63,6 @@ def __init__(self, mix_params): def __str__(self): return self.__class__.__name__+"_"+str(self.mix_params.mass_order) - - def Pmf_HighDensityLimit(self): - """ The probability that a given flavor state is a particular matter state in the - infinite density limit. - - Returns - ------- - 6 x 6 matrix - """ - PmfHDL = np.zeros((6,6)) - - NU_E, NU_MU, NU_TAU, NU_E_BAR, NU_MU_BAR, NU_TAU_BAR = \ - ThreeFlavor.NU_E, ThreeFlavor.NU_MU, ThreeFlavor.NU_TAU, \ - ThreeFlavor.NU_E_BAR, ThreeFlavor.NU_MU_BAR, ThreeFlavor.NU_TAU_BAR - - M2 = np.zeros((6,6)) - M2[1,1] = self.mix_params.dm21_2.value - M2[2,2] = self.mix_params.dm31_2.value - M2[1+3,1+3] = self.mix_params.dm21_2.value - M2[2+3,2+3] = self.mix_params.dm31_2.value - - U = self.mix_params.VacuumMixingMatrix() - - HV = U @ M2 @ np.conjugate(np.transpose(U)) - - T = np.real( ( HV[NU_MU,NU_MU] + HV[NU_TAU,NU_TAU] ) / 2 ) - D = np.abs( HV[NU_MU,NU_TAU] )**2 - np.abs( HV[NU_MU,NU_MU]-HV[NU_TAU,NU_TAU] )**2 / 4 - - Tbar = np.real( ( HV[NU_MU_BAR,NU_MU_BAR] + HV[NU_TAU_BAR,NU_TAU_BAR] ) / 2 ) - Dbar = np.abs( HV[NU_MU_BAR,NU_TAU_BAR] )**2 \ - -np.abs( HV[NU_MU_BAR,NU_MU_BAR]-HV[NU_TAU_BAR,NU_TAU_BAR] )**2 / 4 - - """ The NMO case. Matter state 3 is the electron neutrino, matter state 1bar is the electron - antineutrino. - """ - if self.mix_params.mass_order == MassHierarchy.NORMAL: - k1 = T - np.sqrt(D) - k2 = T + np.sqrt(D) - k2bar = Tbar - np.sqrt(Dbar) - k3bar = Tbar + np.sqrt(Dbar) - - PmfHDL[0,NU_MU] = ( np.real(HV[NU_TAU,NU_TAU]) - k1 ) / ( k2 - k1 ) - PmfHDL[0,NU_TAU] = ( np.real(HV[NU_MU,NU_MU]) - k1 ) / ( k2 - k1 ) - PmfHDL[1,NU_MU] = ( np.real(HV[NU_TAU,NU_TAU]) - k2 ) / ( k1 - k2 ) - PmfHDL[1,NU_TAU] = ( np.real(HV[NU_MU,NU_MU]) - k2 ) / ( k1 - k2 ) - PmfHDL[2,NU_E] = 1 - - PmfHDL[3,NU_E_BAR] = 1 - PmfHDL[4,NU_MU_BAR] = ( np.real(HV[NU_TAU_BAR,NU_TAU_BAR]) - k2bar ) / ( k3bar - k2bar ) - PmfHDL[4,NU_TAU_BAR] = ( np.real(HV[NU_MU_BAR,NU_MU_BAR]) - k2bar ) / ( k3bar - k2bar ) - PmfHDL[5,NU_MU_BAR] = ( np.real(HV[NU_TAU_BAR,NU_TAU_BAR]) - k3bar ) / ( k2bar - k3bar ) - PmfHDL[5,NU_TAU_BAR] = ( np.real(HV[NU_MU_BAR,NU_MU_BAR]) - k3bar ) / ( k2bar - k3bar ) - - """ The IMO case. Matter state 2 is the electron neutrino, matter state 3bar is the electron - antineutrino. - """ - if self.mix_params.mass_order == MassHierarchy.INVERTED: - k1 = T + np.sqrt(D) - k3 = T - np.sqrt(D) - k1bar = Tbar - np.sqrt(Dbar) - k2bar = Tbar + np.sqrt(Dbar) - - PmfHDL[0,NU_MU] = ( np.real(HV[NU_TAU,NU_TAU]) - k1 ) / ( k3 - k1 ) - PmfHDL[0,NU_TAU] = ( np.real(HV[NU_MU,NU_MU]) - k1 ) / ( k3 - k1 ) - PmfHDL[1,NU_E] = 1 - PmfHDL[2,NU_MU] = ( np.real(HV[NU_TAU,NU_TAU]) - k3 ) / ( k1 - k3 ) - PmfHDL[2,NU_TAU] = ( np.real(HV[NU_MU,NU_MU]) - k3 ) / ( k1 - k3 ) - - PmfHDL[3,NU_MU_BAR] = ( np.real(HV[NU_TAU_BAR,NU_TAU_BAR]) - k1bar ) / ( k2bar - k1bar ) - PmfHDL[3,NU_TAU_BAR] = ( np.real(HV[NU_MU_BAR,NU_MU_BAR]) - k1bar ) / ( k2bar - k1bar ) - PmfHDL[4,NU_MU_BAR] = ( np.real(HV[NU_TAU_BAR,NU_TAU_BAR]) - k2bar ) / ( k1bar - k2bar ) - PmfHDL[4,NU_TAU_BAR] = ( np.real(HV[NU_MU_BAR,NU_MU_BAR]) - k2bar ) / ( k1bar - k2bar ) - PmfHDL[5,NU_E_BAR] = 1 - - return PmfHDL ############################################################################### @@ -255,6 +183,7 @@ def get_probabilities(self, t, E): return p def apply(self, flux): + """This transformation returns the object without transform""" return flux ############################################################################### @@ -277,7 +206,7 @@ class AdiabaticMSW(ThreeFlavorTransformation): """Adiabatic MSW effect.""" def get_probabilities(self, t, E): - PHDL = self.Pmf_HighDensityLimit() + PHDL = self.mix_params.Pmf_HighDensityLimit() return FlavorMatrix(PHDL,ThreeFlavor) ############################################################################### diff --git a/python/snewpy/neutrino.py b/python/snewpy/neutrino.py index cc1b6423d..ba222566a 100644 --- a/python/snewpy/neutrino.py +++ b/python/snewpy/neutrino.py @@ -8,6 +8,8 @@ from collections.abc import Mapping from .flavor import ThreeFlavor as Flavor +from .flavor import TwoFlavor, ThreeFlavor +from .flavor import FlavorScheme, FlavorMatrix class MassHierarchy(IntEnum): @@ -53,7 +55,13 @@ class ThreeFlavorMixingParameters(Mapping): #mass ordering mass_order: MassHierarchy | None = None # Note: in IH, the mass splittings are: m3..............m1..m2. - + + #define the basis states + basis_mass = FlavorScheme("ThreeFlavor_MassBasis", start=0, + names=['NU_1','NU_2','NU_3','NU_1_BAR','NU_2_BAR','NU_3_BAR']) + basis_flavor = FlavorScheme("ThreeFlavor_FlavorBasis", start=0, + names=['NU_E','NU_MU','NU_TAU','NU_E_BAR','NU_MU_BAR','NU_TAU_BAR']) + def __iter__(self): return iter(self.__dict__) def __getitem__(self, item): @@ -126,33 +134,83 @@ def VacuumMixingMatrix(self): U = self.ComplexRotationMatrix(1,2,self.theta23,0) \ @ self.ComplexRotationMatrix(0,2,self.theta13,self.deltaCP) \ @ self.ComplexRotationMatrix(0,1,self.theta12,0) - - """Reorder rows to match SNEWPY's flavor ordering convention""" - U[ThreeFlavor.NU_E], U[ThreeFlavor.NU_MU], U[ThreeFlavor.NU_TAU], \ - U[ThreeFlavor.NU_E_BAR], U[ThreeFlavor.NU_MU_BAR], U[ThreeFlavor.NU_TAU_BAR] = \ - U[0], U[1], U[2], U[3], U[4], U[5] - - return U + return FlavorMatrix(U, flavor=self.basis_flavor, from_flavor=self.basis_mass) def ComplexRotationMatrix(self,i,j,theta,phase): """A complex rotation matrix. N.B. the minus sign in the complex exponential matches PDG convention""" - V = np.zeros((6,6),dtype = 'complex_') - for k in range(6): - V[k,k] = 1 - - V[i,i] = np.cos(theta) - V[j,j] = V[i,i] - V[i,j] = np.sin(theta) * ( np.cos(phase) - 1j * np.sin(phase) ) - V[j,i] = - np.conjugate(V[i,j]) - - V[i+3,i+3] = V[i,i] - V[j+3,j+3] = V[j,j] - V[i+3,j+3] = np.conjugate(V[i,j]) - V[j+3,i+3] = np.conjugate(V[j,i]) - + theta = (theta< Date: Thu, 25 Jul 2024 17:58:26 +0300 Subject: [PATCH 18/42] Added _repr_markdown_ for FlavorMatrix --- python/snewpy/flavor.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/python/snewpy/flavor.py b/python/snewpy/flavor.py index 2f56c588a..9484b8173 100644 --- a/python/snewpy/flavor.py +++ b/python/snewpy/flavor.py @@ -114,6 +114,23 @@ def __repr__(self): if(len(self.shape)==2): s+='\n'+repr(self.array) return s + + def _repr_markdown_(self): + if(len(self.shape)>2): + return self.__repr__() + #make a markdown table + res = [f'**{self.__class__.__name__}**:<`{self.flavor_in.__name__}`->`{self.flavor_out.__name__}`> shape={self.shape}',''] + flavors0 = [f.to_tex() for f in self.flavor_in] + flavors1 = [f.to_tex() for f in self.flavor_out] + hdr = '|'.join(['*']+flavors1) + res+=[hdr] + res+=['|'.join(['-:',]*(len(flavors1)+1))] + for f0 in self.flavor_in: + line = [f0.to_tex()]+[f'{v:.3g}' for v in self[:,f0]] + res+=['|'.join(line)] + res+=['---------'] + return '\n'.join(res) + def __eq__(self,other): return self.flavor_in==other.flavor_in and self.flavor_out==other.flavor_out and np.allclose(self.array,other.array) From f04d49510e227a9d6f7ca220d6d6770fdacc2a04 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Thu, 25 Jul 2024 18:09:31 +0300 Subject: [PATCH 19/42] Corrected ComplexRotationMatrix creation --- python/snewpy/neutrino.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/python/snewpy/neutrino.py b/python/snewpy/neutrino.py index ba222566a..7662306ef 100644 --- a/python/snewpy/neutrino.py +++ b/python/snewpy/neutrino.py @@ -141,7 +141,7 @@ def ComplexRotationMatrix(self,i,j,theta,phase): """A complex rotation matrix. N.B. the minus sign in the complex exponential matches PDG convention""" theta = (theta< Date: Thu, 25 Jul 2024 23:41:12 +0300 Subject: [PATCH 20/42] Added '__mul__','real','imag','abs' and 'abs2' methods for FlavorMatrix --- python/snewpy/flavor.py | 21 +++++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/python/snewpy/flavor.py b/python/snewpy/flavor.py index 9484b8173..b2a6b76ce 100644 --- a/python/snewpy/flavor.py +++ b/python/snewpy/flavor.py @@ -133,7 +133,24 @@ def _repr_markdown_(self): def __eq__(self,other): return self.flavor_in==other.flavor_in and self.flavor_out==other.flavor_out and np.allclose(self.array,other.array) - + @property + def real(self): + return FlavorMatrix(self.array.real, self.flavor_out, from_flavor = self.flavor_in) + @property + def imag(self): + return FlavorMatrix(self.array.imag, self.flavor_out, from_flavor = self.flavor_in) + + def abs(self): + return FlavorMatrix(np.abs(self.array), self.flavor_out, from_flavor = self.flavor_in) + def abs2(self): + return FlavorMatrix(np.abs(self.array**2), self.flavor_out, from_flavor = self.flavor_in) + + def __mul__(self, other): + if isinstance(other, FlavorMatrix): + if not ((other.flavor_in==self.flavor_in)and(other.flavor_out==self.flavor_out)): + raise TypeError(f"Cannot multiply matrices with different flavor schemes: {self._repr_short()} and {other._repr_short()}") + other = other.array + return FlavorMatrix(self.array*other, self.flavor_out, from_flavor = self.flavor_in) def __matmul__(self, other): if isinstance(other, FlavorMatrix): assert self.flavor_in==other.flavor_out, f"Incompatible spaces {self.flavor_in}!={other.flavor_out}" @@ -174,7 +191,7 @@ def conjugate(self): "apply complex conjugate" return FlavorMatrix(array = self.array.conjugate(), flavor = self.flavor_out, from_flavor=self.flavor_in) - + @classmethod def eye(cls, flavor:FlavorScheme, from_flavor:FlavorScheme = None): from_flavor = from_flavor or flavor From 4a604f033d6a3cc41d0c654482928087df8dfea9 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Fri, 26 Jul 2024 13:14:32 +0300 Subject: [PATCH 21/42] Implement AdiabaticMSW, NonAdiabaticMSWH, ThreeFlavorDecoherence, TwoFlavorDecoherence --- python/snewpy/flavor_transformation.py | 163 ++++++++----------------- 1 file changed, 50 insertions(+), 113 deletions(-) diff --git a/python/snewpy/flavor_transformation.py b/python/snewpy/flavor_transformation.py index d3a62d0ea..4b94803ac 100644 --- a/python/snewpy/flavor_transformation.py +++ b/python/snewpy/flavor_transformation.py @@ -52,7 +52,7 @@ def P(self, t, E): class ThreeFlavorTransformation(FlavorTransformation): """Base class defining common data and methods for all three flavor transformations""" - def __init__(self, mix_params): + def __init__(self, mix_params=None): """Initialize flavor transformation Parameters @@ -205,95 +205,58 @@ def get_probabilities(self, t, E): class AdiabaticMSW(ThreeFlavorTransformation): """Adiabatic MSW effect.""" - def get_probabilities(self, t, E): - PHDL = self.mix_params.Pmf_HighDensityLimit() - return FlavorMatrix(PHDL,ThreeFlavor) + def get_probabilities(self, t, E): + Pmf = self.mix_params.Pmf_HighDensityLimit() + + D = self.mix_params.VacuumMixingMatrix().abs2() + return D@Pmf + ############################################################################### class NonAdiabaticMSWH(ThreeFlavorTransformation): - """Nonadiabatic MSW H resonance""" - - def _get_SNprobabilities(self, t, E): - """neutrino and antineutrino transition probabilities in the supernova - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - Pmf : 6 x 6 array - """ - PHDL = self.Pmf_HighDensityLimit() - Pmf = np.empty((r,6,len(E))) - - for m in range(len(E)): - Pmf[:,:,m] = PHDL[:,:] - if self.mix_params.mass_order == MassHierarchy.NORMAL: - for f in ThreeFlavor: - Pmf[1,f,m], Pmf[2,f,m] = Pmf[2,f,m], Pmf[1,f,m] - if self.mix_params.mass_order == MassHierarchy.INVERTED: - for f in ThreeFlavor: - Pmf[3,f,m], Pmf[5,f,m] = Pmf[5,f,m], Pmf[3,f,m] - - return Pmf - - def get_probabilities(self, t, E): - Pmf = self._get_SNprobabilities(t,E) - D = ThreeFlavorNoEarthMatter(self.mix_params).get_probabilities(t,E) - - p = np.empty((6,6,len(E))) - for m in range(len(E)): - p[:,:,m] = D[:,:,m] @ Pmf[:,:,m] + """Nonadiabatic MSW H resonance. + The NonAdiabaticMSWH transformation assumes that the H resonance mixing is nonadiabatic. + This case is relevant when a shock is present at the H resonance densities (Schirato & Fuller 2002). + + For the NMO the H resonance occurs in the neutrinos (Kneller & McLaughlin 2009) between ‘matter’ + states ν2 and ν3. + In the IMO the H resonance mixes the antineutrino matter states ν̄1 and ν̄3. + """ + def get_probabilities(self, t, E): + Pmf = self.mix_params.Pmf_HighDensityLimit() + if self.mix_params.mass_order == MassHierarchy.NORMAL: + Pmf[['NU_2','NU_3'],:] = Pmf[['NU_3','NU_2'],:] + else: + Pmf[['NU_1_BAR','NU_3_BAR'],:] = Pmf[['NU_3_BAR','NU_1_BAR'],:] - return p + D = self.mix_params.VacuumMixingMatrix().abs2() + return D@Pmf ############################################################################### class TwoFlavorDecoherence(ThreeFlavorTransformation): - """equal mixing of whatever two matter states form the MSW H resonance""" - - def get_SNprobabilities(self, t, E): - """neutrino and antineutrino transition probabilities. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - p : 6 x 6 array or array of 6 x 6 arrays - """ - PHDL = self.Pmf_HighDensityLimit(E) - Pmf = np.empty((6,6,len(E))) - - for m in range(len(E)): - Pmf[:,:,m] = PHDL[:,:] - if self.mix_params.mass_order == MassHierarchy.NORMAL: - for f in ThreeFlavor: - Pmf[1,f,m], Pmf[2,f,m] = 0.5 * ( Pmf[1,f,m] + Pmf[2,f,m] ), 0.5 * ( Pmf[1,f,m] + Pmf[2,f,m] ) - if self.mix_params.mass_order == MassHierarchy.INVERTED: - for f in ThreeFlavor: - Pmf[3,f,m], Pmf[5,f,m] = 0.5 * ( Pmf[3,f,m] + Pmf[5,f,m] ), 0.5 * ( Pmf[3,f,m] + Pmf[5,f,m] ) - - return Pmf + """equal mixing of whatever two matter states form the MSW H resonance. + The TwoFlavorDecoherence transformation is relevant when the size of the density fluctuations + is ≲ 10% for densities around the H resonance density —see Kneller (2010); Kneller & Mauney (2013). + This prescription models the effect of the turbulence as leading to 50% mixing between the + matter states which participate in the H resonance. + + In the NMO this is ν2 and ν3, + For the IMO, the H resonance occurs in the antineutrinos between antineutrino matter states ν̄1 +and ν̄3. + """ + def get_probabilities(self, t, E): - Pmf = self.get_SNprobabilities(t,E) - D = ThreeFlavorNoEarthMatter(self.mix_params).get_probabilities(t,E) - - p = np.empty((6,6,len(E))) - for m in range(len(E)): - p[:,:,m] = D[:,:,m] @ Pmf[:,:,m] - - return p + Pmf = self.mix_params.Pmf_HighDensityLimit() + if self.mix_params.mass_order == MassHierarchy.NORMAL: + Pmf['NU_2']=Pmf['NU_3']=0.5*(Pmf['NU_2']+Pmf['NU_3']) + else: + Pmf['NU_1_BAR']=Pmf['NU_3_BAR']=0.5*(Pmf['NU_1_BAR']+Pmf['NU_3_BAR']) + + D = self.mix_params.VacuumMixingMatrix().abs2() + return D@Pmf ############################################################################### @@ -307,27 +270,12 @@ def __init__(self): def __str__(self): return f'ThreeFlavorDecoherence' - def get_SNprobabilities(self, t, E): - """neutrino and antineutrino transition probabilities. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - p : 6 x 6 array or array of 6 x 6 arrays - """ - p = np.ones((6,6,len(E))) / 3 - - return p - def get_probabilities(self, t, E): """Equal mixing so Earth matter has no effect""" - return self.get_SNprobabilities(t,E) + @FlavorMatrix.from_function(ThreeFlavor) + def Prob_ff(f1,f2): + return (f1.is_neutrino==f2.is_neutrino)*1/3. + return Prob_ff ############################################################################### @@ -361,20 +309,9 @@ def __init__(self, SNprofile, mix_params = None, rmin = None, rmax = None): rmax : the ending radius of the calculation. rmax will be corrected by SNOSHEWS to the maximum radius of the profile if rmax is greater than that value """ self.SNprofile = SNprofile - - super().__init__(mix_params) - - if rmin == None: - rmin = 0 - else: - self.rmin = rmin - - if rmax == None: - rmax = 1e99 - else: - self.rmax = rmax - - super().__init__(mix_params) + super().__init__(mix_params) + self.rmin = rmin or 0 + self.rmax = rmax or 1e99 def __str__(self): return f'MSW_' + str(self.mix_params.mass_order) @@ -432,7 +369,7 @@ def get_SNprobabilities(self, t, E): pSN = SNOSHEWS.Run(ID) # restructure the results - Pmf = np.zeros((6,6,ID.NE)) + Pmf = FlavorMatrix(np.zeros((6,6,ID.NE)) for m in range(ID.NE): Pmf[0,ThreeFlavor.NU_E,m] = pSN[0][m][0][0] From 7620e78f622d1e0434a51f76591f6a94ac2c6568 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Fri, 26 Jul 2024 13:40:43 +0300 Subject: [PATCH 22/42] Converting the thansformation matrix to the flux flavor scheme --- python/snewpy/flavor_transformation.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/python/snewpy/flavor_transformation.py b/python/snewpy/flavor_transformation.py index 4b94803ac..6b60b5a5e 100644 --- a/python/snewpy/flavor_transformation.py +++ b/python/snewpy/flavor_transformation.py @@ -43,6 +43,7 @@ def get_probabilities(self, t, E)->FlavorMatrix: def apply(self, flux): M = self.get_probabilities(flux.time, flux.energy) + M = (flux.flavor_scheme< Date: Fri, 26 Jul 2024 13:41:19 +0300 Subject: [PATCH 23/42] Correct check for the case 'flavor_out is None' --- python/snewpy/flavor.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/python/snewpy/flavor.py b/python/snewpy/flavor.py index b2a6b76ce..c184a6107 100644 --- a/python/snewpy/flavor.py +++ b/python/snewpy/flavor.py @@ -207,7 +207,8 @@ def zeros(cls, flavor:FlavorScheme, from_flavor:FlavorScheme = None): @classmethod def from_function(cls, flavor:FlavorScheme, from_flavor:FlavorScheme = None): """A decorator for creating the flavor matrix from the given function""" - from_flavor = from_flavor or flavor + if from_flavor is None: + from_flavor = flavor def _decorator(function): data = [[function(f1,f2) for f2 in from_flavor] From c0cb00b963b402dddb21c5bdcb7349b9ea567816 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Sat, 27 Jul 2024 20:25:20 +0300 Subject: [PATCH 24/42] Fix the CompleteExchange matrix --- python/snewpy/flavor_transformation.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/python/snewpy/flavor_transformation.py b/python/snewpy/flavor_transformation.py index 6b60b5a5e..f58e079dc 100644 --- a/python/snewpy/flavor_transformation.py +++ b/python/snewpy/flavor_transformation.py @@ -194,12 +194,11 @@ class CompleteExchange(FlavorTransformation): are half exchanged with the mu flavors and the half with the tau flavors. """ def get_probabilities(self, t, E): - p = FlavorMatrix.zeros(TwoFlavor) - p['NU_X','NU_E'] = 1 - p['NU_E','NU_X'] = 1 - p['NU_E_BAR','NU_X_BAR'] = 1 - p['NU_X_BAR','NU_E_BAR'] = 1 - return (ThreeFlavor< Date: Sat, 27 Jul 2024 23:34:25 +0300 Subject: [PATCH 25/42] Adding FourFlavorMixingParameters.Pmf_HDL --- python/snewpy/neutrino.py | 107 ++++++++++++++++++++++++++++++-------- 1 file changed, 86 insertions(+), 21 deletions(-) diff --git a/python/snewpy/neutrino.py b/python/snewpy/neutrino.py index 7662306ef..29642d68e 100644 --- a/python/snewpy/neutrino.py +++ b/python/snewpy/neutrino.py @@ -234,7 +234,12 @@ class FourFlavorMixingParameters(ThreeFlavorMixingParameters): dm41_2: u.Quantity[u.eV**2] = 0< Date: Sun, 28 Jul 2024 00:00:52 +0300 Subject: [PATCH 26/42] Implemented AdiabaticMSWes/NonAdiabaticMSWes --- python/snewpy/flavor_transformation.py | 247 ++----------------------- 1 file changed, 17 insertions(+), 230 deletions(-) diff --git a/python/snewpy/flavor_transformation.py b/python/snewpy/flavor_transformation.py index f58e079dc..c04574d35 100644 --- a/python/snewpy/flavor_transformation.py +++ b/python/snewpy/flavor_transformation.py @@ -67,113 +67,18 @@ def __str__(self): ############################################################################### -class FourFlavorTransformation: +class FourFlavorTransformation(ThreeFlavorTransformation): """Base class defining common data and method for all four flavor transformations""" - def __init__(self, mix_params): + def __init__(self, mix_params=None): """Initialize flavor transformation Parameters ---------- mix_params : FourFlavorMixingParameters instance """ - if mix_params == None: - self.mix_params = FourFlavorMixingParameters(MassHierarchy.NORMAL) - else: - self.mix_params = mix_params - - def __str__(self): - return self.__class__.__name__+"_"+str(self.mix_params.mass_order) - - def Pmf_HighDensityLimit(self): - """ The probability that a given flavor state is 'detected' as a particular matter state in the - infinite density limit. This method assumes the 4th matter state is the - heaviest and that the electron fraction remains larger than 1/3 - - Returns - ------- - 8 x 8 matrix - """ - PmfHDL = np.zeros((8,8)) - - NU_E, NU_MU, NU_TAU, NU_S, NU_E_BAR, NU_MU_BAR, NU_TAU_BAR, NU_S_BAR = \ - FourFlavor.NU_E, FourFlavor.NU_MU, FourFlavor.NU_TAU, FourFlavor.NU_S, \ - FourFlavor.NU_E_BAR, FourFlavor.NU_MU_BAR, FourFlavor.NU_TAU_BAR, FourFlavor.NU_S_BAR - - M2 = np.zeros((8,8)) - M2[1,1] = self.mix_params.dm21_2.value - M2[2,2] = self.mix_params.dm31_2.value - M2[3,3] = self.mix_params.dm41_2.value - M2[1+4,1+4] = self.mix_params.dm21_2.value - M2[2+4,2+4] = self.mix_params.dm31_2.value - M2[3+4,3+4] = self.mix_params.dm41_2.value - - U = VacuumMixingMatrix() - - HV = np.zeros((6,6,len(E))) - for m in range(len(E)): - HV[:,:,m] = U @ M2 @ np.conjugate(np.transpose(U)) / ( 2 * E[m] ) - - T = np.real( (HV[NU_MU,NU_MU] + HV[NU_TAU,NU_TAU]) / 2 ) - D = np.abs( HV[NU_MU,NU_TAU] )**2 - np.abs( HV[NU_MU,NU_MU]-HV[NU_TAU,NU_TAU] )**2 / 4 - - Tbar = np.real( ( HV[NU_MU_BAR,NU_MU_BAR] + HV[NU_TAU_BAR,NU_TAU_BAR] ) / 2 ) - Dbar = np.abs( HV[NU_MU_BAR,NU_TAU_BAR] )**2 \ - -np.abs( HV[NU_MU_BAR,NU_MU_BAR]-HV[NU_TAU_BAR,NU_TAU_BAR] )**2 / 4 - - U = Vacuum_Mixing_Matrix() - - """ The NMO case. Matter state 4 is the electron neutrino, matter state 1bar is the electron - antineutrino. Matter state 3 is the sterile neutrino, matter state 2bar is the sterile - antineutrino. - """ - if self.mix_params.mass_order == MassHierarchy.NORMAL: - k1 = T - np.sqrt(D) - k2 = T + np.sqrt(D) - k3bar = Tbar - np.sqrt(Dbar) - k4bar = Tbar + np.sqrt(Dbar) - - PmfHDL[0,NU_MU] = ( np.real(HV[NU_TAU,NU_TAU]) - k1 ) / ( k2 - k1 ) - PmfHDL[0,NU_TAU] = ( np.real(HV[NU_MU,NU_MU]) - k1 ) / ( k2 - k1 ) - PmfHDL[1,NU_MU] = ( np.real(HV[NU_TAU,NU_TAU]) - k2 ) / ( k1 - k2 ) - PmfHDL[1,NU_TAU] = ( np.real(HV[NU_MU,NU_MU]) - k2 ) / ( k1 - k2 ) - PmfHDL[2,NU_S] = 1 - PmfHDL[3,NU_E] = 1 - - PmfHDL[4,NU_E_BAR] = 1 - PmfHDL[5,NU_S_BAR] = 1 - PmfHDL[6,NU_MU_BAR] = ( np.real(HV[NU_TAU_BAR,NU_TAU_BAR]) - k3bar ) / ( k4bar - k3bar ) - PmfHDL[6,NU_TAU_BAR] = ( np.real(HV[NU_MU_BAR,NU_MU_BAR]) - k3bar ) / ( k4bar - k3bar ) - PmfHDL[7,NU_MU_BAR] = ( np.real(HV[NU_TAU_BAR,NU_TAU_BAR]) - k4bar ) / ( k3bar - k4bar ) - PmfHDL[7,NU_TAU_BAR] = ( np.real(HV[NU_MU_BAR,NU_MU_BAR]) - k4bar ) / ( k3bar - k4bar ) - - - """ The IMO case. Matter state 4 is the electron neutrino, matter state 3bar is the electron - antineutrino. Matter state 2 is the sterile neutrino, matter state 1bar is the sterile - antineutrino. - """ - if self.mix_params.mass_order == MassHierarchy.INVERTED: - k1 = T + np.sqrt(D) - k3 = T - np.sqrt(D) - k2bar = Tbar - np.sqrt(Dbar) - k4bar = Tbar + np.sqrt(Dbar) - - PmfHDL[0,NU_MU] = ( np.real(HV[NU_TAU,NU_TAU]) - k1 ) / ( k3 - k1 ) - PmfHDL[0,NU_TAU] = ( np.real(HV[NU_MU,NU_MU]) - k1 ) / ( k3 - k1 ) - PmfHDL[1,NU_S] = 1 - PmfHDL[2,NU_MU] = ( np.real(HV[NU_TAU,NU_TAU]) - k3 ) / ( k1 - k3 ) - PmfHDL[2,NU_TAU] = ( np.real(HV[NU_MU,NU_MU]) - k3 ) / ( k1 - k3 ) - PmfHDL[3,NU_E] = 1 - - PmfHDL[4,NU_S_BAR] = 1 - PmfHDL[5,NU_MU_BAR] = ( np.real(HV[NU_TAU_BAR,NU_TAU_BAR]) - k2bar ) / ( k4bar - k2bar ) - PmfHDL[5,NU_TAU_BAR] = ( np.real(HV[NU_MU_BAR,NU_MU_BAR]) - k2bar ) / ( k4bar - k2bar ) - PmfHDL[6,NU_E_BAR] = 1 - PmfHDL[7,NU_TAU_BAR] = ( np.real(HV[NU_MU_BAR,NU_MU_BAR]) - k4bar ) / ( k2bar - k4bar ) - PmfHDL[7,NU_MU_BAR] = ( np.real(HV[NU_TAU_BAR,NU_TAU_BAR]) - k4bar ) / ( k2bar - k4bar ) - - return PmfHDL - + self.mix_params = mix_params or FourFlavorMixingParameters(**MixingParameters(MassHierarchy.NORMAL)) + ############################################################################### class NoTransformation(FlavorTransformation): @@ -433,69 +338,10 @@ class AdiabaticMSWes(FourFlavorTransformation): For further insight see, for example, Esmaili, Peres, and Serpico, Phys. Rev. D 90, 033013 (2014). """ - def __init__(self, mix_params = None): - """Initialize flavor transformation - - Parameters - ---------- - mix_params : FourFlavorMixingParameters instance or None - """ - super().__init__(mix_params) - - def __str__(self): - return f'AdiabaticMSWes_' + str(self.mix_params.mass_order) - - def get_SNprobabilities(self, t, E): - """neutrino and antineutrino transition probabilities. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - Pmf : 8 x 8 matrix - """ - PHDL = self.Pmf_HighDensityLimit() - Pmf = np.empty((8,8,len(E))) - - for m in range(len(E)): - Pmf[:,:,m] = PHDL[:,:] - - return Pmf - - - def get_probabilities(self, t, E): - """neutrino and antineutrino transition probabilities. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - p : 6 x 6 matrix - """ - p = np.empty((8,8,len(E))) - - D = FourFlavorNoEarthMatter(self.mix_params).get_probabilities(t,E) - Pmf = self.get_SNprobabilities(t,E) - - # multiply the D matrix and the Pmf matrix together - for m in range(len(E)): - p[:,:,m] = D[:,:,m] @ Pmf[:,:,m] - - # remove sterile rows/columns: A is a 6 x 8 matrix - A = Remove_Steriles() - p = A @ p @ np.transpose(A) - - return p + def get_probabilities(self, t, E): + Pmf = self.mix_params.Pmf_HighDensityLimit() + D = self.mix_params.VacuumMixingMatrix().abs2() + return D@Pmf ############################################################################### @@ -511,77 +357,18 @@ class NonAdiabaticMSWes(FourFlavorTransformation): For further insight see, for example, Esmaili, Peres, and Serpico, Phys. Rev. D 90, 033013 (2014). """ - def __init__(self, mix_params = None): - """Initialize flavor transformation - - Parameters - ---------- - mix_params : FourFlavorMixingParameters instance or None - """ - super().__init__(mix_params) - - def __str__(self): - return f'NonAdiabaticMSWes_' + str(self.mix_params.mass_order) - - def get_SNprobabilities(self, t, E): - """neutrino and antineutrino transition probabilities. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - Pmf : an 8 x 8 matrix - """ - PHDL = self.Pmf_HighDensityLimit() - Pmf = np.empty((8,8,len(E))) - - for m in range(len(E)): - Pmf[:,:,m] = PHDL[:,:] - - if self.mix_params.mass_order == MassHierarchy.NORMAL: - for f in ThreeFlavor: - Pmf[2,f,m], Pmf[3,f,m] = Pmf[3,f,m], Pmf[2,f,m] - Pmf[5,f,m], Pmf[6,f,m], Pmf[7,f,m] = Pmf[6,f,m], Pmf[7,f,m], Pmf[5,f,m] - - if self.mix_params.mass_order == MassHierarchy.INVERTED: - for f in ThreeFlavor: - Pmf[1,f,m], Pmf[3,f,m] = Pmf[3,f,m], Pmf[1,f,m] - Pmf[4,f,m], Pmf[5,f,m], Pmf[7,f,m] = Pmf[5,f,m], Pmf[7,f,m], Pmf[4,f,m] - - return Pmf def get_probabilities(self, t, E): - """neutrino and antineutrino transition probabilities. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - p : 6 x 6 array or array of 6 x 6 arrays - """ - - D = FourFlavorNoEarthMatter(self.mix_params).get_probabilities(t,E) - Pmf = self.get_SNprobabilities(t,E) - - p = np.empty((8,8,len(E))) - for m in range(len(E)): - p[:,:,m] = D[:,:,m] @ Pmf[:,:,m] + Pmf = self.mix_params.Pmf_HighDensityLimit() + if self.mix_params.mass_order == MassHierarchy.NORMAL: + Pmf[['NU_3','NU_4'],:] = Pmf[['NU_4','NU_3'],:] + Pmf[['NU_2_BAR','NU_3_BAR','NU_4_BAR'],:] = Pmf[['NU_3_BAR','NU_4_BAR','NU_2_BAR'],:] + else: + Pmf[['NU_2','NU_4'],:] = Pmf[['NU_4','NU_2'],:] + Pmf[['NU_1_BAR','NU_2_BAR','NU_4_BAR'],:] = Pmf[['NU_2_BAR','NU_4_BAR','NU_1_BAR'],:] - # remove sterile rows/columns - A = Remove_Steriles() - p = A @ p @ np.transpose(A) - - return p + D = self.mix_params.VacuumMixingMatrix().abs2() + return D@Pmf ############################################################################### From de73410f65479a6ef773309f44781c1d5a0d2317 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Mon, 29 Jul 2024 11:05:29 +0300 Subject: [PATCH 27/42] Adding 'project_to' method --- python/snewpy/flux.py | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/python/snewpy/flux.py b/python/snewpy/flux.py index 511c1e296..82c2ba9e6 100644 --- a/python/snewpy/flux.py +++ b/python/snewpy/flux.py @@ -189,7 +189,7 @@ def __getitem__(self, args)->'Container': if isinstance(args[0],str) or isinstance(args[0],FlavorScheme): args[0] = self.flavor_scheme[args[0]] for n,arg in enumerate(args): - if not isinstance(arg, slice): + if isinstance(arg, int): arg = slice(arg, arg + 1) arg_slices[n] = arg @@ -445,12 +445,21 @@ def __class_getitem__(cls, args): cls._unit_classes[unit] = type(name,(cls,),{'unit':unit}) return cls._unit_classes[unit] - def plot(flux, projection='energy', styles=None, **kwargs): - if projection=='energy': - fP = flux.integrate_or_sum('time') + def project_to(self, axis='energy', squeeze=False): + if axis=='energy': + fP = self.integrate_or_sum('time') + else: + fP = self.integrate_or_sum('energy') + x = fP.__dict__[axis] + if squeeze: + return x, fP.array.squeeze().T else: - fP = flux.integrate_or_sum('energy') - x = fP.__dict__[projection] + return x, fP + + + def plot(flux, projection='energy', styles=None, **kwargs): + x, fP = flux.project_to(projection, squeeze=False) + if isinstance(styles, dict): styles = styles.get elif styles==None: From ad0aecbde076a22c6776be96363810c6428f24af Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Mon, 29 Jul 2024 15:28:40 +0300 Subject: [PATCH 28/42] Fix the presn test --- python/snewpy/test/test_presn_rates.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/python/snewpy/test/test_presn_rates.py b/python/snewpy/test/test_presn_rates.py index 719e81dac..8f82a25ce 100644 --- a/python/snewpy/test/test_presn_rates.py +++ b/python/snewpy/test/test_presn_rates.py @@ -1,6 +1,7 @@ import astropy.units as u import numpy as np from snewpy.models import presn +from snewpy.neutrino import MixingParameters from snewpy.flavor_transformation import AdiabaticMSW, MassHierarchy from snewpy.rate_calculator import RateCalculator import pytest @@ -16,7 +17,7 @@ @pytest.mark.xfail(reason="`model.get_flux` uses `get_transformed_flux`, which is currently hard coded to a TwoFlavor scheme.", raises=AttributeError) @pytest.mark.parametrize('model_class',[presn.Odrzywolek_2010, presn.Kato_2017, presn.Patton_2017, presn.Yoshida_2016]) -@pytest.mark.parametrize('transformation',[AdiabaticMSW(mh=mh) for mh in MassHierarchy]) +@pytest.mark.parametrize('transformation',[AdiabaticMSW(MixingParameters(mh)) for mh in MassHierarchy]) @pytest.mark.parametrize('detector', ["wc100kt30prct"]) def test_presn_rate(model_class, transformation, detector): model = model_class(progenitor_mass=15*u.Msun) From ab35a44697518183d59628ed8bc8e63e14ce5108 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Mon, 29 Jul 2024 15:30:49 +0300 Subject: [PATCH 29/42] Allowing a shornt notation in FlavorScheme - without 'NU_' prefix --- python/snewpy/flavor.py | 5 +++++ python/snewpy/test/test_flavors.py | 14 ++++++++++++++ 2 files changed, 19 insertions(+) diff --git a/python/snewpy/flavor.py b/python/snewpy/flavor.py index c184a6107..2381bd13b 100644 --- a/python/snewpy/flavor.py +++ b/python/snewpy/flavor.py @@ -17,6 +17,11 @@ def __getitem__(cls, key): return key #if this is a string find it by name if isinstance(key, str): + #prepare the string + key=key.upper() + if not key.startswith('NU_'): + key = 'NU_'+key + #try to get the proper values try: return super().__getitem__(key) except KeyError as e: diff --git a/python/snewpy/test/test_flavors.py b/python/snewpy/test/test_flavors.py index ead9e0dc6..cb55f5e04 100644 --- a/python/snewpy/test/test_flavors.py +++ b/python/snewpy/test/test_flavors.py @@ -17,8 +17,14 @@ def test_flavor_scheme_lengths(): def test_getitem_string(): assert TwoFlavor['NU_E'] == TwoFlavor.NU_E assert TwoFlavor['NU_X'] == TwoFlavor.NU_X + #short notations + assert ThreeFlavor['E'] == ThreeFlavor['e'] == ThreeFlavor.NU_E + assert ThreeFlavor['MU'] == ThreeFlavor['mu'] == ThreeFlavor.NU_MU + assert ThreeFlavor['MU_BAR'] == ThreeFlavor['mu_bar'] == ThreeFlavor.NU_MU_BAR with pytest.raises(KeyError): TwoFlavor['NU_MU'] + with pytest.raises(KeyError): + ThreeFlavor['NU_X'] @staticmethod def test_getitem_enum(): @@ -100,6 +106,14 @@ def test_getitem(): assert np.allclose(m['NU_E'], m['NU_E',:]) assert np.allclose(m[:,:], m.array) + @staticmethod + def test_getitem_short(): + m = FlavorMatrix.eye(ThreeFlavor,ThreeFlavor) + assert m['NU_E','NU_E']==m['e','e'] + assert m['NU_MU','NU_E']==m['mu','e'] + assert m['NU_E_BAR','NU_E']==m['e_bar','e'] + assert m['NU_TAU_BAR','NU_TAU']==m['tau_bar','tau'] + @staticmethod def test_setitem(): m = FlavorMatrix.eye(TwoFlavor,TwoFlavor) From e40c0477058394faa855b98429b7f5d61daa79c3 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Mon, 29 Jul 2024 19:26:37 +0300 Subject: [PATCH 30/42] Using shorter flavor notations in formulas --- python/snewpy/neutrino.py | 116 +++++++++++++++++++------------------- 1 file changed, 57 insertions(+), 59 deletions(-) diff --git a/python/snewpy/neutrino.py b/python/snewpy/neutrino.py index 29642d68e..68dd9bc48 100644 --- a/python/snewpy/neutrino.py +++ b/python/snewpy/neutrino.py @@ -165,12 +165,12 @@ def Pmf_HighDensityLimit(self): U = self.VacuumMixingMatrix() HV = U @ M2 @ np.conjugate(U.T) - T = np.real( ( HV['NU_MU','NU_MU'] + HV['NU_TAU','NU_TAU'] ) / 2 ) - D = np.abs( HV['NU_MU','NU_TAU'] )**2 \ - -np.abs( HV['NU_MU','NU_MU']-HV['NU_TAU','NU_TAU'] )**2 / 4 - Tbar = np.real( ( HV["NU_MU_BAR","NU_MU_BAR"] + HV["NU_TAU_BAR","NU_TAU_BAR"] ) / 2 ) - Dbar = np.abs( HV["NU_MU_BAR","NU_TAU_BAR"] )**2 \ - -np.abs( HV["NU_MU_BAR","NU_MU_BAR"]-HV["NU_TAU_BAR","NU_TAU_BAR"] )**2 / 4 + T = np.real( ( HV['mu','mu'] + HV['tau','tau'] ) / 2 ) + D = np.abs( HV['mu','tau'] )**2 \ + -np.abs( HV['mu','mu']-HV['tau','tau'] )**2 / 4 + Tbar = np.real( ( HV['mu_bar','mu_bar'] + HV['tau_bar','tau_bar'] ) / 2 ) + Dbar = np.abs( HV['mu_bar','tau_bar'] )**2 \ + -np.abs( HV['mu_bar','mu_bar']-HV['tau_bar','tau_bar'] )**2 / 4 PmfHDL = FlavorMatrix.zeros(self.basis_mass, self.basis_flavor) @@ -180,17 +180,17 @@ def Pmf_HighDensityLimit(self): k2 = T + np.sqrt(D) k2bar = Tbar - np.sqrt(Dbar) k3bar = Tbar + np.sqrt(Dbar) - PmfHDL['NU_1','NU_MU'] = (HV['NU_TAU','NU_TAU'].real - k1)/(k2-k1) - PmfHDL['NU_1','NU_TAU'] = (HV['NU_MU','NU_MU'].real - k1)/(k2-k1) - PmfHDL['NU_2','NU_MU'] = (HV['NU_TAU','NU_TAU'].real - k2)/(k1-k2) - PmfHDL['NU_2','NU_TAU'] = (HV['NU_MU','NU_MU'].real - k2)/(k1 -k2) - PmfHDL['NU_3','NU_E'] = 1 + PmfHDL['1','mu'] = (HV['tau','tau'].real - k1)/(k2-k1) + PmfHDL['1','tau'] = (HV['mu','mu'].real - k1)/(k2-k1) + PmfHDL['2','mu'] = (HV['tau','tau'].real - k2)/(k1-k2) + PmfHDL['2','tau'] = (HV['mu','mu'].real - k2)/(k1 -k2) + PmfHDL['3','e'] = 1 - PmfHDL['NU_1_BAR','NU_E_BAR'] = 1 - PmfHDL['NU_2_BAR','NU_MU_BAR'] = (HV['NU_TAU_BAR','NU_TAU_BAR'].real - k2bar ) / ( k3bar - k2bar ) - PmfHDL['NU_2_BAR','NU_TAU_BAR'] = (HV['NU_MU_BAR','NU_MU_BAR'].real - k2bar ) / ( k3bar - k2bar ) - PmfHDL['NU_3_BAR','NU_MU_BAR'] = (HV['NU_TAU_BAR','NU_TAU_BAR'].real - k3bar ) / ( k2bar - k3bar ) - PmfHDL['NU_3_BAR','NU_TAU_BAR'] = (HV['NU_MU_BAR','NU_MU_BAR'].real - k3bar ) / ( k2bar - k3bar ) + PmfHDL['1_bar','e_bar'] = 1 + PmfHDL['2_bar','mu_bar'] = (HV['tau_bar','tau_bar'].real - k2bar ) / ( k3bar - k2bar ) + PmfHDL['2_bar','tau_bar'] = (HV['mu_bar','mu_bar'].real - k2bar ) / ( k3bar - k2bar ) + PmfHDL['3_bar','mu_bar'] = (HV['tau_bar','tau_bar'].real - k3bar ) / ( k2bar - k3bar ) + PmfHDL['3_bar','tau_bar'] = (HV['mu_bar','mu_bar'].real - k3bar ) / ( k2bar - k3bar ) elif self.mass_order == MassHierarchy.INVERTED: # The IMO case. Matter state 2 is the electron neutrino, matter state 3bar is the electron antineutrino. @@ -199,17 +199,17 @@ def Pmf_HighDensityLimit(self): k1bar = Tbar - np.sqrt(Dbar) k2bar = Tbar + np.sqrt(Dbar) - PmfHDL['NU_1','NU_MU'] = ( HV['NU_TAU','NU_TAU'].real - k1 ) / ( k3 - k1 ) - PmfHDL['NU_1','NU_TAU'] = ( HV['NU_MU','NU_MU'].real - k1 ) / ( k3 - k1 ) - PmfHDL['NU_2','NU_E'] = 1 - PmfHDL['NU_3','NU_MU'] = ( HV['NU_TAU','NU_TAU'].real - k3) / ( k1 - k3 ) - PmfHDL['NU_3','NU_TAU'] = ( HV['NU_MU','NU_MU'].real - k3 ) / ( k1 - k3 ) + PmfHDL['1','mu'] = ( HV['tau','tau'].real - k1 ) / ( k3 - k1 ) + PmfHDL['1','tau'] = ( HV['mu','mu'].real - k1 ) / ( k3 - k1 ) + PmfHDL['2','e'] = 1 + PmfHDL['3','mu'] = ( HV['tau','tau'].real - k3) / ( k1 - k3 ) + PmfHDL['3','tau'] = ( HV['mu','mu'].real - k3 ) / ( k1 - k3 ) - PmfHDL['NU_1_BAR','NU_MU_BAR'] = ( HV['NU_TAU_BAR','NU_TAU_BAR'].real - k1bar ) / ( k2bar - k1bar ) - PmfHDL['NU_1_BAR','NU_TAU_BAR'] = ( HV['NU_MU_BAR','NU_MU_BAR'].real - k1bar ) / ( k2bar - k1bar ) - PmfHDL['NU_2_BAR','NU_MU_BAR'] = ( HV['NU_TAU_BAR','NU_TAU_BAR'].real - k2bar ) / ( k1bar - k2bar ) - PmfHDL['NU_2_BAR','NU_TAU_BAR'] = ( HV['NU_MU_BAR','NU_MU_BAR'].real - k2bar ) / ( k1bar - k2bar ) - PmfHDL['NU_3_BAR','NU_E_BAR'] = 1 + PmfHDL['1_bar','mu_bar'] = ( HV['tau_bar','tau_bar'].real - k1bar ) / ( k2bar - k1bar ) + PmfHDL['1_bar','tau_bar'] = ( HV['mu_bar','mu_bar'].real - k1bar ) / ( k2bar - k1bar ) + PmfHDL['2_bar','mu_bar'] = ( HV['tau_bar','tau_bar'].real - k2bar ) / ( k1bar - k2bar ) + PmfHDL['2_bar','tau_bar'] = ( HV['mu_bar','mu_bar'].real - k2bar ) / ( k1bar - k2bar ) + PmfHDL['3_bar','e_bar'] = 1 return PmfHDL @dataclass @@ -236,9 +236,9 @@ class FourFlavorMixingParameters(ThreeFlavorMixingParameters): dm43_2: u.Quantity | None = None #define the basis states - basis_mass = FlavorScheme("ThreeFlavor_MassBasis", start=0, + basis_mass = FlavorScheme("FourFlavor_MassBasis", start=0, names=['NU_1','NU_2','NU_3','NU_4','NU_1_BAR','NU_2_BAR','NU_3_BAR','NU_4_BAR']) - basis_flavor = FlavorScheme("ThreeFlavor_FlavorBasis", start=0, + basis_flavor = FlavorScheme("FourFlavor_FlavorBasis", start=0, names=['NU_E','NU_MU','NU_TAU','NU_S','NU_E_BAR','NU_MU_BAR','NU_TAU_BAR','NU_S_BAR']) def __post_init__(self): super().__post_init__() @@ -315,13 +315,11 @@ def Pmf_HighDensityLimit(self): U = self.VacuumMixingMatrix() HV = U @ M2 @ np.conjugate(U.T) - T = np.real( ( HV['NU_MU','NU_MU'] + HV['NU_TAU','NU_TAU'] ) / 2 ) - D = np.abs( HV['NU_MU','NU_TAU'] )**2 \ - -np.abs( HV['NU_MU','NU_MU']-HV['NU_TAU','NU_TAU'] )**2 / 4 - Tbar = np.real( ( HV["NU_MU_BAR","NU_MU_BAR"] + HV["NU_TAU_BAR","NU_TAU_BAR"] ) / 2 ) - Dbar = np.abs( HV["NU_MU_BAR","NU_TAU_BAR"] )**2 \ - -np.abs( HV["NU_MU_BAR","NU_MU_BAR"]-HV["NU_TAU_BAR","NU_TAU_BAR"] )**2 / 4 - + T = np.real( ( HV['mu','mu'] + HV['tau','tau'] ) / 2 ) + D = np.abs(HV['mu','tau'])**2-np.abs( HV['mu','mu']-HV['tau','tau'] )**2 / 4 + Tbar = np.real( ( HV['mu_bar','mu_bar'] + HV['tau_bar','tau_bar'] ) / 2 ) + Dbar = np.abs( HV['mu_bar','tau_bar'] )**2 \ + -np.abs( HV['mu_bar','mu_bar']-HV['tau_bar','tau_bar'] )**2 / 4 PmfHDL = FlavorMatrix.zeros(self.basis_mass, self.basis_flavor) if self.mass_order == MassHierarchy.NORMAL: @@ -332,19 +330,19 @@ def Pmf_HighDensityLimit(self): k3bar = Tbar - np.sqrt(Dbar) k4bar = Tbar + np.sqrt(Dbar) - PmfHDL['NU_1','NU_MU'] = (HV['NU_TAU','NU_TAU'].real - k1)/(k2-k1) - PmfHDL['NU_1','NU_TAU'] = (HV['NU_MU','NU_MU'].real - k1)/(k2-k1) - PmfHDL['NU_2','NU_MU'] = (HV['NU_TAU','NU_TAU'].real - k2)/(k1-k2) - PmfHDL['NU_2','NU_TAU'] = (HV['NU_MU','NU_MU'].real - k2)/(k1-k2) - PmfHDL['NU_3','NU_S'] = 1 - PmfHDL['NU_4','NU_E'] = 1 + PmfHDL['1','mu'] = (HV['tau','tau'].real - k1)/(k2-k1) + PmfHDL['1','tau'] = (HV['mu','mu'].real - k1)/(k2-k1) + PmfHDL['2','mu'] = (HV['tau','tau'].real - k2)/(k1-k2) + PmfHDL['2','tau'] = (HV['mu','mu'].real - k2)/(k1-k2) + PmfHDL['3','S'] = 1 + PmfHDL['4','e'] = 1 - PmfHDL['NU_1_BAR','NU_E_BAR'] = 1 - PmfHDL['NU_2_BAR','NU_S_BAR'] = 1 - PmfHDL['NU_3_BAR','NU_MU_BAR'] = (HV['NU_TAU_BAR','NU_TAU_BAR'].real - k3bar ) / ( k4bar - k3bar ) - PmfHDL['NU_3_BAR','NU_TAU_BAR'] = (HV['NU_MU_BAR','NU_MU_BAR'].real - k3bar ) / ( k4bar - k3bar ) - PmfHDL['NU_4_BAR','NU_MU_BAR'] = (HV['NU_TAU_BAR','NU_TAU_BAR'].real - k4bar ) / ( k3bar - k4bar ) - PmfHDL['NU_4_BAR','NU_TAU_BAR'] = (HV['NU_MU_BAR','NU_MU_BAR'].real - k4bar ) / ( k3bar - k4bar ) + PmfHDL['1_bar','e_bar'] = 1 + PmfHDL['2_bar','S_bar'] = 1 + PmfHDL['3_bar','mu_bar'] = (HV['tau_bar','tau_bar'].real - k3bar ) / ( k4bar - k3bar ) + PmfHDL['3_bar','tau_bar'] = (HV['mu_bar','mu_bar'].real - k3bar ) / ( k4bar - k3bar ) + PmfHDL['4_bar','mu_bar'] = (HV['tau_bar','tau_bar'].real - k4bar ) / ( k3bar - k4bar ) + PmfHDL['4_bar','tau_bar'] = (HV['mu_bar','mu_bar'].real - k4bar ) / ( k3bar - k4bar ) elif self.mass_order == MassHierarchy.INVERTED: # The IMO case. Matter state 4 is the electron neutrino, matter state 3bar is the electron antineutrino. @@ -354,19 +352,19 @@ def Pmf_HighDensityLimit(self): k2bar = Tbar - np.sqrt(Dbar) k4bar = Tbar + np.sqrt(Dbar) - PmfHDL['NU_1','NU_MU'] = ( HV['NU_TAU','NU_TAU'].real - k1 ) / ( k3 - k1 ) - PmfHDL['NU_1','NU_TAU'] = ( HV['NU_MU','NU_MU'].real - k1 ) / ( k3 - k1 ) - PmfHDL['NU_2','NU_S'] = 1 - PmfHDL['NU_3','NU_MU'] = ( HV['NU_TAU','NU_TAU'].real - k3) / ( k1 - k3 ) - PmfHDL['NU_3','NU_TAU'] = ( HV['NU_MU','NU_MU'].real - k3 ) / ( k1 - k3 ) - PmfHDL['NU_4','NU_E'] = 1 + PmfHDL['1','mu'] = ( HV['tau','tau'].real - k1 ) / ( k3 - k1 ) + PmfHDL['1','tau'] = ( HV['mu','mu'].real - k1 ) / ( k3 - k1 ) + PmfHDL['2','S'] = 1 + PmfHDL['3','mu'] = ( HV['tau','tau'].real - k3) / ( k1 - k3 ) + PmfHDL['3','tau'] = ( HV['mu','mu'].real - k3 ) / ( k1 - k3 ) + PmfHDL['4','e'] = 1 - PmfHDL['NU_1_BAR','NU_S_BAR'] = 1 - PmfHDL['NU_2_BAR','NU_MU_BAR'] = ( HV['NU_TAU_BAR','NU_TAU_BAR'].real - k2bar ) / ( k4bar - k2bar ) - PmfHDL['NU_2_BAR','NU_TAU_BAR'] = ( HV['NU_MU_BAR','NU_MU_BAR'].real - k2bar ) / ( k4bar - k2bar ) - PmfHDL['NU_3_BAR','NU_E_BAR'] = 1 - PmfHDL['NU_4_BAR','NU_MU_BAR'] = ( HV['NU_TAU_BAR','NU_TAU_BAR'].real - k4bar ) / ( k2bar - k4bar ) - PmfHDL['NU_4_BAR','NU_TAU_BAR'] = ( HV['NU_MU_BAR','NU_MU_BAR'].real - k4bar ) / ( k2bar - k4bar ) + PmfHDL['1_bar','S_bar'] = 1 + PmfHDL['2_bar','mu_bar'] = ( HV['tau_bar','tau_bar'].real - k2bar ) / ( k4bar - k2bar ) + PmfHDL['2_bar','tau_bar'] = ( HV['mu_bar','mu_bar'].real - k2bar ) / ( k4bar - k2bar ) + PmfHDL['3_bar','e_bar'] = 1 + PmfHDL['4_bar','mu_bar'] = ( HV['tau_bar','tau_bar'].real - k4bar ) / ( k2bar - k4bar ) + PmfHDL['4_bar','tau_bar'] = ( HV['mu_bar','mu_bar'].real - k4bar ) / ( k2bar - k4bar ) return PmfHDL parameter_presets = { From 1380b633f81ec12add5e9cc3c31753dcc3880af9 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Mon, 29 Jul 2024 22:46:50 +0300 Subject: [PATCH 31/42] Update flavor_transformations tests (partially) --- python/snewpy/test/test_04_xforms.py | 678 +++++++++++++-------------- 1 file changed, 317 insertions(+), 361 deletions(-) diff --git a/python/snewpy/test/test_04_xforms.py b/python/snewpy/test/test_04_xforms.py index 0ed0f99e3..595252612 100644 --- a/python/snewpy/test/test_04_xforms.py +++ b/python/snewpy/test/test_04_xforms.py @@ -1,380 +1,336 @@ # -*- coding: utf-8 -*- import unittest +import pytest -from snewpy.flavor_transformation import MassHierarchy, MixingParameters -from snewpy.flavor_transformation import \ - NoTransformation, CompleteExchange, \ - AdiabaticMSW, NonAdiabaticMSWH, \ - AdiabaticMSWes, NonAdiabaticMSWes, \ - TwoFlavorDecoherence, ThreeFlavorDecoherence, \ - NeutrinoDecay, QuantumDecoherence +from snewpy.neutrino import MassHierarchy, MixingParameters, FourFlavorMixingParameters +import snewpy.flavor_transformation as xforms + +from snewpy.flavor import ThreeFlavor, FourFlavor, TwoFlavor from astropy import units as u from astropy import constants as c import numpy as np -from numpy import sin, cos, exp, abs - -class TestFlavorTransformations(unittest.TestCase): - - def setUp(self): - # Dummy time and energy arrays, with proper dimensions. - self.t = np.arange(10) * u.s - self.E = np.linspace(1,100,21) * u.MeV - - # Dummy mixing angles. - self.theta12 = 33 * u.deg - self.theta13 = 9 * u.deg - self.theta23 = 49 * u.deg - self.theta14 = 1 * u.deg - - # Dummy neutrino decay parameters; see arXiv:1910.01127. - self.mass3 = 0.5 * u.eV/c.c**2 - self.lifetime = 1 * u.day - self.distance = 10 * u.kpc - - # Quantum decoherence parameters; see arXiv:2306.17591. - self.gamma3 = (1e-27 * u.eV / (c.hbar * c.c)).to('1/kpc') - self.gamma8 = (1e-27 * u.eV / (c.hbar * c.c)).to('1/kpc') - self.n = 0 - self.energy_ref = 10 * u.MeV - - def test_noxform(self): +from numpy import sin,cos,abs + +# Dummy neutrino decay parameters; see arXiv:1910.01127. +@pytest.fixture +def params_NeutrinoDecay(): + return { + 'mass3':0.5 * u.eV/c.c**2, + 'lifetime':1 * u.day, + 'distance':10 * u.kpc + } + +@pytest.fixture +def params_QuantumDecoherence(): + return { + 'gamma3': (1e-27 * u.eV / (c.hbar * c.c)).to('1/kpc'), + 'gamma8': (1e-27 * u.eV / (c.hbar * c.c)).to('1/kpc'), + 'n': 0, + 'energy_ref': 10 * u.MeV + } + +@pytest.fixture(autouse=True) +def t(): + return np.arange(10) * u.s +@pytest.fixture(autouse=True) +def E(): + return np.linspace(1,100,21) * u.MeV +class TestFlavorTransformations: + + def test_NoTransformation(self): """ Survival probabilities for no oscillations """ - xform = NoTransformation() - - self.assertEqual(xform.prob_ee(self.t, self.E), 1) - self.assertEqual(xform.prob_ex(self.t, self.E), 0) - self.assertEqual(xform.prob_xx(self.t, self.E), 1) - self.assertEqual(xform.prob_xe(self.t, self.E), 0) - - self.assertEqual(xform.prob_eebar(self.t, self.E), 1) - self.assertEqual(xform.prob_exbar(self.t, self.E), 0) - self.assertEqual(xform.prob_xxbar(self.t, self.E), 1) - self.assertEqual(xform.prob_xebar(self.t, self.E), 0) - - def test_fullex(self): + P = xforms.NoTransformation().P(t, E) + assert P.flavor_in == ThreeFlavor + assert P.flavor_out == ThreeFlavor + for f1 in ThreeFlavor: + for f2 in ThreeFlavor: + assert np.isclose(P[f1,f2] , 0 if f1!=f2 else 1) + + def test_CompleteExchange(self): """ Survival probabilities for complete electron->X transformation """ - xform = CompleteExchange() + P = xforms.CompleteExchange().P(t, E) + assert P.flavor_in == ThreeFlavor + assert P.flavor_out == ThreeFlavor + #general check + for f1 in ThreeFlavor: + for f2 in ThreeFlavor: + if f1.is_neutrino!=f2.is_neutrino: + assert np.isclose(P[f1,f2],0) + else: + assert np.isclose(P[f1,f2],0.5*(f1!=f2)) - self.assertEqual(xform.prob_ee(self.t, self.E), 0) - self.assertEqual(xform.prob_ex(self.t, self.E), 1) - self.assertEqual(xform.prob_xx(self.t, self.E), 0.5) - self.assertEqual(xform.prob_xe(self.t, self.E), 0.5) - - self.assertEqual(xform.prob_eebar(self.t, self.E), 0) - self.assertEqual(xform.prob_exbar(self.t, self.E), 1) - self.assertEqual(xform.prob_xxbar(self.t, self.E), 0.5) - self.assertEqual(xform.prob_xebar(self.t, self.E), 0.5) - - def test_adiabaticmsw_nmo(self): + def test_AdiabaticMSW_NMO(self): """ Adiabatic MSW with normal ordering """ - xform = AdiabaticMSW(mix_angles=(self.theta12, self.theta13, self.theta23), mh=MassHierarchy.NORMAL) - - self.assertEqual(xform.prob_ee(self.t, self.E), sin(self.theta13)**2) - self.assertEqual(xform.prob_ex(self.t, self.E), 1. - sin(self.theta13)**2) - self.assertEqual(xform.prob_xx(self.t, self.E), 0.5*(1. + sin(self.theta13)**2)) - self.assertEqual(xform.prob_xe(self.t, self.E), 0.5*(1. - sin(self.theta13)**2)) - - self.assertEqual(xform.prob_eebar(self.t, self.E), (cos(self.theta12)*cos(self.theta13))**2) - self.assertEqual(xform.prob_exbar(self.t, self.E), 1. - (cos(self.theta12)*cos(self.theta13))**2) - self.assertEqual(xform.prob_xxbar(self.t, self.E), 0.5*(1. + (cos(self.theta12)*cos(self.theta13))**2)) - self.assertEqual(xform.prob_xebar(self.t, self.E), 0.5*(1. - (cos(self.theta12)*cos(self.theta13))**2)) - - # Test interface using default mixing angles defined in the submodule. mixpars = MixingParameters(MassHierarchy.NORMAL) th12, th13, th23 = mixpars.get_mixing_angles() + + P = xforms.AdiabaticMSW(mixpars).P(t, E) + #convert to TwoFlavor case + P = (TwoFlavor<10% turbulence and are now disabled - # in the class. - self.assertFalse(xform.prob_ee(self.t, self.E) == 0.5) - self.assertFalse(xform.prob_ex(self.t, self.E) == 0.5) - self.assertFalse(xform.prob_xx(self.t, self.E) == 0.75) - self.assertFalse(xform.prob_xe(self.t, self.E) == 0.25) - - self.assertFalse(xform.prob_eebar(self.t, self.E) == 0.5) - self.assertFalse(xform.prob_exbar(self.t, self.E) == 0.5) - self.assertFalse(xform.prob_xxbar(self.t, self.E) == 0.75) - self.assertFalse(xform.prob_xebar(self.t, self.E) == 0.25) - # Calculation that applies to 1% turbulence. - mixpars = MixingParameters() + mixpars = MixingParameters('NORMAL') th12, th13, th23 = mixpars.get_mixing_angles() + + P = xforms.TwoFlavorDecoherence(mixpars).P(t,E) + #convert to TwoFlavor case + P = (TwoFlavor< Date: Mon, 29 Jul 2024 23:03:19 +0300 Subject: [PATCH 32/42] Using cdouble data type instead of complex_ --- python/snewpy/neutrino.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/snewpy/neutrino.py b/python/snewpy/neutrino.py index 68dd9bc48..5d79fab29 100644 --- a/python/snewpy/neutrino.py +++ b/python/snewpy/neutrino.py @@ -141,7 +141,7 @@ def ComplexRotationMatrix(self,i,j,theta,phase): """A complex rotation matrix. N.B. the minus sign in the complex exponential matches PDG convention""" theta = (theta< Date: Mon, 29 Jul 2024 23:13:35 +0300 Subject: [PATCH 33/42] Remove erroneous 'bar_bar' indices --- python/snewpy/test/test_04_xforms.py | 46 ++++++++++++++-------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/python/snewpy/test/test_04_xforms.py b/python/snewpy/test/test_04_xforms.py index 595252612..634f0d9e6 100644 --- a/python/snewpy/test/test_04_xforms.py +++ b/python/snewpy/test/test_04_xforms.py @@ -328,9 +328,9 @@ def test_ThreeFlavorDecoherence(self): assert (abs(P['x','e'] - 1./3) < 1e-12) assert np.isclose(P['e_bar','e_bar'] , 1./3) - assert (abs(P['e_bar','xbar_bar'] - 2./3) < 1e-12) - assert (abs(P['x_bar','xbar_bar'] - 2./3) < 1e-12) - assert (abs(P['x_bar','ebar_bar'] - 1./3) < 1e-12) + assert (abs(P['e_bar','x_bar'] - 2./3) < 1e-12) + assert (abs(P['x_bar','x_bar'] - 2./3) < 1e-12) + assert (abs(P['x_bar','e_bar'] - 1./3) < 1e-12) def test_nudecay_nmo(self): """ @@ -358,8 +358,8 @@ def test_nudecay_nmo(self): prob_exbar = np.asarray([De1*(1.-exp(-xform.gamma(_E)*self.distance)) + De2 + De3*exp(-xform.gamma(_E)*self.distance) for _E in self.E]) assert np.isclose(P['e_bar','e_bar'], De3) - assert (np.array_equal(P['e_bar','xbar_bar'], prob_exbar)) - assert (np.array_equal(P['x_bar','xbar_bar'], 1. - 0.5*prob_exbar)) + assert (np.array_equal(P['e_bar','x_bar'], prob_exbar)) + assert (np.array_equal(P['x_bar','x_bar'], 1. - 0.5*prob_exbar)) assert np.isclose(P['x_bar','e_bar'], 0.5*(1. - De3)) def test_nudecay_nmo_default_mixing(self): @@ -387,8 +387,8 @@ def test_nudecay_nmo_default_mixing(self): prob_exbar = np.asarray([De1*(1.-exp(-xform.gamma(_E)*self.distance)) + De2 + De3*exp(-xform.gamma(_E)*self.distance) for _E in self.E]) assert np.isclose(P['e_bar','e_bar'], De3) - assert (np.array_equal(P['e_bar','xbar_bar'], prob_exbar)) - assert (np.array_equal(P['x_bar','xbar_bar'], 1. - 0.5*prob_exbar)) + assert (np.array_equal(P['e_bar','x_bar'], prob_exbar)) + assert (np.array_equal(P['x_bar','x_bar'], 1. - 0.5*prob_exbar)) assert np.isclose(P['x_bar','e_bar'], 0.5*(1. - De3)) def test_nudecay_imo(self): @@ -414,8 +414,8 @@ def test_nudecay_imo(self): De3*(1-np.exp(-xform.gamma(_E)*self.distance)) for _E in self.E]) assert np.isclose(P['e_bar','e_bar'], De3) - assert (np.array_equal(P['e_bar','xbar_bar'], prob_exbar)) - assert (np.array_equal(P['x_bar','xbar_bar'], 1. - 0.5*prob_exbar)) + assert (np.array_equal(P['e_bar','x_bar'], prob_exbar)) + assert (np.array_equal(P['x_bar','x_bar'], 1. - 0.5*prob_exbar)) assert np.isclose(P['x_bar','e_bar'], 0.5*(1. - De3)) def test_nudecay_imo_default_mixing(self): @@ -445,8 +445,8 @@ def test_nudecay_imo_default_mixing(self): De3*(1-np.exp(-xform.gamma(_E)*self.distance)) for _E in self.E]) assert np.isclose(P['e_bar','e_bar'], De3) - assert (np.array_equal(P['e_bar','xbar_bar'], prob_exbar)) - assert (np.array_equal(P['x_bar','xbar_bar'], 1. - 0.5*prob_exbar)) + assert (np.array_equal(P['e_bar','x_bar'], prob_exbar)) + assert (np.array_equal(P['x_bar','x_bar'], 1. - 0.5*prob_exbar)) assert np.isclose(P['x_bar','e_bar'], 0.5*(1. - De3)) def test_quantumdecoherence_nmo(self): @@ -479,9 +479,9 @@ def test_quantumdecoherence_nmo(self): prob_eebar = np.asarray([xform.P11(_E)*De1 + xform.P21(_E)*De2 + xform.P31(_E)*De3 for _E in self.E]) - assert (np.array_equal(P['e_bar','xbar_bar'], 1 - prob_eebar)) - assert (np.array_equal(P['x_bar','xbar_bar'], 1. - 0.5*(1 - prob_eebar))) - assert (np.array_equal(P['x_bar','ebar_bar'], 0.5*(1. - prob_eebar))) + assert (np.array_equal(P['e_bar','x_bar'], 1 - prob_eebar)) + assert (np.array_equal(P['x_bar','x_bar'], 1. - 0.5*(1 - prob_eebar))) + assert (np.array_equal(P['x_bar','e_bar'], 0.5*(1. - prob_eebar))) def test_quantumdecoherence_nmo_default_mixing(self): """ @@ -507,9 +507,9 @@ def test_quantumdecoherence_nmo_default_mixing(self): prob_eebar = np.asarray([xform.P11(_E)*De1 + xform.P21(_E)*De2 + xform.P31(_E)*De3 for _E in self.E]) - assert (np.array_equal(P['e_bar','xbar_bar'], 1 - prob_eebar)) - assert (np.array_equal(P['x_bar','xbar_bar'], 1. - 0.5*(1 - prob_eebar))) - assert (np.array_equal(P['x_bar','ebar_bar'], 0.5*(1. - prob_eebar))) + assert (np.array_equal(P['e_bar','x_bar'], 1 - prob_eebar)) + assert (np.array_equal(P['x_bar','x_bar'], 1. - 0.5*(1 - prob_eebar))) + assert (np.array_equal(P['x_bar','e_bar'], 0.5*(1. - prob_eebar))) def test_quantumdecoherence_imo(self): """ @@ -541,9 +541,9 @@ def test_quantumdecoherence_imo(self): prob_eebar = np.asarray([xform.P31(_E)*De1 + xform.P32(_E)*De2 + xform.P33(_E)*De3 for _E in self.E]) - assert (np.array_equal(P['e_bar','xbar_bar'], 1 - prob_eebar)) - assert (np.array_equal(P['x_bar','xbar_bar'], 1. - 0.5*(1 - prob_eebar))) - assert (np.array_equal(P['x_bar','ebar_bar'], 0.5*(1. - prob_eebar))) + assert (np.array_equal(P['e_bar','x_bar'], 1 - prob_eebar)) + assert (np.array_equal(P['x_bar','x_bar'], 1. - 0.5*(1 - prob_eebar))) + assert (np.array_equal(P['x_bar','e_bar'], 0.5*(1. - prob_eebar))) def test_quantumdecoherence_imo_default_mixing(self): """ @@ -569,6 +569,6 @@ def test_quantumdecoherence_imo_default_mixing(self): prob_eebar = np.asarray([xform.P31(_E)*De1 + xform.P32(_E)*De2 + xform.P33(_E)*De3 for _E in self.E]) - assert (np.array_equal(P['e_bar','xbar_bar'], 1 - prob_eebar)) - assert (np.array_equal(P['x_bar','xbar_bar'], 1. - 0.5*(1 - prob_eebar))) - assert (np.array_equal(P['x_bar','ebar_bar'], 0.5*(1. - prob_eebar))) \ No newline at end of file + assert (np.array_equal(P['e_bar','x_bar'], 1 - prob_eebar)) + assert (np.array_equal(P['x_bar','x_bar'], 1. - 0.5*(1 - prob_eebar))) + assert (np.array_equal(P['x_bar','e_bar'], 0.5*(1. - prob_eebar))) \ No newline at end of file From 9852dbd7ece066ba60ebcf0cd591bf7d8cbac0a5 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Tue, 30 Jul 2024 17:24:40 +0300 Subject: [PATCH 34/42] Adding (flavor1<FlavorMatrix: + if isinstance(obj, EnumMeta): + return conversion_matrix(from_flavor=flv,to_flavor=obj) + elif hasattr(obj, '__lshift__'): + return obj<FlavorMatrix: + if isinstance(obj, EnumMeta): + return conversion_matrix(from_flavor=obj,to_flavor=flv) + elif hasattr(obj, '__rshift__'): + return obj>>flv + else: + raise TypeError(f'Cannot apply flavor conversion to object of type {type(obj)}') + + FlavorScheme.conversion_matrix = classmethod(conversion_matrix) -EnumMeta.__rshift__ = conversion_matrix -EnumMeta.__lshift__ = lambda f1,f2:conversion_matrix(f2,f1) +EnumMeta.__rshift__ = rshift +EnumMeta.__lshift__ = lshift diff --git a/python/snewpy/test/test_flavors.py b/python/snewpy/test/test_flavors.py index cb55f5e04..28fc26b37 100644 --- a/python/snewpy/test/test_flavors.py +++ b/python/snewpy/test/test_flavors.py @@ -141,7 +141,7 @@ def test_conversion_matrices_for_same_flavor_are_unity(): matrix = flavor>>flavor assert isinstance(matrix, FlavorMatrix) assert np.allclose(matrix.array, np.eye(len(flavor))) - + @staticmethod @pytest.mark.parametrize('flavor_in',flavor_schemes) @pytest.mark.parametrize('flavor_out',flavor_schemes) @@ -150,4 +150,24 @@ def test_conversion_matrices(flavor_in, flavor_out): assert M==flavor_out< Date: Tue, 30 Jul 2024 18:42:10 +0300 Subject: [PATCH 35/42] Adding 'extra_dims;' argument to FlavorMatrix.zeros and FlavorMatrix.eye to create larger matrices --- python/snewpy/flavor.py | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/python/snewpy/flavor.py b/python/snewpy/flavor.py index 2f0de6cce..a82a31439 100644 --- a/python/snewpy/flavor.py +++ b/python/snewpy/flavor.py @@ -198,15 +198,12 @@ def conjugate(self): flavor = self.flavor_out, from_flavor=self.flavor_in) @classmethod - def eye(cls, flavor:FlavorScheme, from_flavor:FlavorScheme = None): - from_flavor = from_flavor or flavor - shape = (len(from_flavor), len(flavor)) - data = np.eye(*shape) - return cls(data, flavor, from_flavor) + def eye(cls, flavor:FlavorScheme, from_flavor:FlavorScheme = None, extra_dims=[]): + return cls.from_function(flavor,from_flavor)(lambda f1,f2: f1==f2*np.ones(shape=extra_dims)) @classmethod - def zeros(cls, flavor:FlavorScheme, from_flavor:FlavorScheme = None): + def zeros(cls, flavor:FlavorScheme, from_flavor:FlavorScheme = None, extra_dims=[]): from_flavor = from_flavor or flavor - shape = (len(from_flavor), len(flavor)) + shape = (len(from_flavor), len(flavor), *extra_dims) data = np.zeros(shape) return cls(data, flavor, from_flavor) @classmethod From 2727d0040fb50dfc812995401a5317c2248d091c Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Tue, 30 Jul 2024 18:54:43 +0300 Subject: [PATCH 36/42] Implement NeutrinoDecay --- python/snewpy/flavor_transformation.py | 29 +++++++++----------------- 1 file changed, 10 insertions(+), 19 deletions(-) diff --git a/python/snewpy/flavor_transformation.py b/python/snewpy/flavor_transformation.py index c04574d35..f947d1ee9 100644 --- a/python/snewpy/flavor_transformation.py +++ b/python/snewpy/flavor_transformation.py @@ -400,9 +400,6 @@ def __init__(self, mix_params=None, mass=1*u.eV/c.c**2, tau=1*u.day, dist=10*u.k self.tau = tau self.d = dist - def __str__(self): - return f'NeutrinoDecay_' + str(self.mix_params.mass_order) - def gamma(self, E): """Decay width of the heaviest neutrino mass. @@ -434,25 +431,19 @@ def get_probabilities(self, t, E): ------- PND : an array of 6 x 6 matrices """ - decay_factor = np.exp(-self.gamma(E)*self.d) - PND = np.zeros((6,6,len(E))) - + decay_factor = np.exp(-self.gamma(E)*self.d) + + PND = FlavorMatrix.eye(self.mix_params.basis_mass, extra_dims=[len(E)]) + if self.mix_params.mass_order == MassHierarchy.NORMAL: - PND[0,0] = 1 - PND[0,2] = 1 - decay_factor - PND[1,1] = 1 - PND[2,2] = decay_factor + PND['NU_1','NU_3'] = 1 - decay_factor + PND['NU_3','NU_3'] = decay_factor - if self.mix_params.mass_order == MassHierarchy.INVERTED: - PND[0,0] = 1 - PND[1,1] = decay_factor - PND[2,1] = 1 - decay_factor - PND[2,2] = 1 - - for i in range(3): - for j in range(3): - PND[i+3,j+3] = PND[i,j] + else: + PND['NU_2','NU_2'] = decay_factor + PND['NU_3','NU_2'] = 1 - decay_factor + PND[3:,3:] = PND[:3,:3] return PND ############################################################################### From dd2c6a872268109c8b95814e01d0e7912597cf6f Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Tue, 30 Jul 2024 19:21:07 +0300 Subject: [PATCH 37/42] Using << operators in tests --- python/snewpy/test/test_04_xforms.py | 112 ++++++++++++++------------- 1 file changed, 59 insertions(+), 53 deletions(-) diff --git a/python/snewpy/test/test_04_xforms.py b/python/snewpy/test/test_04_xforms.py index 634f0d9e6..62cd8274e 100644 --- a/python/snewpy/test/test_04_xforms.py +++ b/python/snewpy/test/test_04_xforms.py @@ -10,7 +10,7 @@ from astropy import units as u from astropy import constants as c import numpy as np -from numpy import sin,cos,abs +from numpy import sin,cos,exp,abs # Dummy neutrino decay parameters; see arXiv:1910.01127. @pytest.fixture @@ -32,12 +32,14 @@ def params_QuantumDecoherence(): @pytest.fixture(autouse=True) def t(): - return np.arange(10) * u.s + return np.arange(10) << u.s + @pytest.fixture(autouse=True) def E(): - return np.linspace(1,100,21) * u.MeV -class TestFlavorTransformations: + return np.linspace(1,100,21) << u.MeV + +class TestFlavorTransformations: def test_NoTransformation(self): """ Survival probabilities for no oscillations @@ -73,7 +75,7 @@ def test_AdiabaticMSW_NMO(self): P = xforms.AdiabaticMSW(mixpars).P(t, E) #convert to TwoFlavor case - P = (TwoFlavor< Date: Wed, 31 Jul 2024 12:22:23 +0300 Subject: [PATCH 38/42] Allow FlavorMatrix slicing return another FlavorMatrix, of the first two dimensions are : --- python/snewpy/flavor.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/python/snewpy/flavor.py b/python/snewpy/flavor.py index a82a31439..8de21987d 100644 --- a/python/snewpy/flavor.py +++ b/python/snewpy/flavor.py @@ -102,11 +102,20 @@ def __init__(self, def _convert_index(self, index): if isinstance(index, str) or (not isinstance(index,typing.Iterable)): index = [index] + #convert flavor dimensions new_idx = [flavors[idx] for idx,flavors in zip(index, [self.flavor_out, self.flavor_in])] + #add remaining dimensions + new_idx+=list(index[2:]) return tuple(new_idx) def __getitem__(self, index): - return self.array[self._convert_index(index)] + index = self._convert_index(index) + array = self.array[index] + #if requested all the flavors, then return a FlavorMatrix using the remaining index + if index[:2]==(slice(None),slice(None)): + return FlavorMatrix(array, self.flavor_out, self.flavor_in) + else: + return array def __setitem__(self, index, value): self.array[self._convert_index(index)] = value From dc02045ea43e629e4db40d86fb27dbe7f0ec41b2 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Wed, 31 Jul 2024 13:36:18 +0300 Subject: [PATCH 39/42] Add check if the conversion matrix is zero --- python/snewpy/flavor.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/python/snewpy/flavor.py b/python/snewpy/flavor.py index 8de21987d..9de15d7cb 100644 --- a/python/snewpy/flavor.py +++ b/python/snewpy/flavor.py @@ -255,6 +255,9 @@ def convert(f1, f2): # convert from more flavors to TwoFlavor return 0.5 return 0. + #check if the conversion matrix is not zero + if np.allclose(convert.array, 0): + raise RuntimeError(f"Conversion matrix {convert._repr_short()} is zero!") return convert def rshift(flv:FlavorScheme, obj:FlavorScheme|FlavorMatrix)->FlavorMatrix: From 5871f51a42d51958cd5ac656c141b7cba0171de7 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Wed, 31 Jul 2024 15:15:58 +0300 Subject: [PATCH 40/42] Adding 'NU' and 'NU_BAR' collective aliases for the FlavorSchemes --- python/snewpy/flavor.py | 12 ++++++++++++ python/snewpy/test/test_flavors.py | 16 +++++++++++++--- 2 files changed, 25 insertions(+), 3 deletions(-) diff --git a/python/snewpy/flavor.py b/python/snewpy/flavor.py index 9de15d7cb..137053ecf 100644 --- a/python/snewpy/flavor.py +++ b/python/snewpy/flavor.py @@ -19,6 +19,12 @@ def __getitem__(cls, key): if isinstance(key, str): #prepare the string key=key.upper() + #check cases for neutrinos and antineutrinos + if key=='NU': + return tuple([f for f in cls.__members__.values() if f.is_neutrino]) + elif key=='NU_BAR': + return tuple([f for f in cls.__members__.values() if not f.is_neutrino]) + #add the prefix if needed if not key.startswith('NU_'): key = 'NU_'+key #try to get the proper values @@ -104,6 +110,12 @@ def _convert_index(self, index): index = [index] #convert flavor dimensions new_idx = [flavors[idx] for idx,flavors in zip(index, [self.flavor_out, self.flavor_in])] + #try to convert first index to list of bracketed index + try: + new_idx[0] = [[f] for f in new_idx[0]] + except: + #it was not iterable + pass #add remaining dimensions new_idx+=list(index[2:]) return tuple(new_idx) diff --git a/python/snewpy/test/test_flavors.py b/python/snewpy/test/test_flavors.py index 28fc26b37..78d30e62c 100644 --- a/python/snewpy/test/test_flavors.py +++ b/python/snewpy/test/test_flavors.py @@ -25,7 +25,12 @@ def test_getitem_string(): TwoFlavor['NU_MU'] with pytest.raises(KeyError): ThreeFlavor['NU_X'] - + @staticmethod + def test_getitem_collective_names(): + assert ThreeFlavor['NU']==(ThreeFlavor.NU_E, ThreeFlavor.NU_MU, ThreeFlavor.NU_TAU) + assert ThreeFlavor['NU']==ThreeFlavor['e','mu','tau'] + assert ThreeFlavor['NU_BAR']==(ThreeFlavor.NU_E_BAR, ThreeFlavor.NU_MU_BAR, ThreeFlavor.NU_TAU_BAR) + assert ThreeFlavor['NU_BAR']==ThreeFlavor['e_bar','mu_bar','tau_bar'] @staticmethod def test_getitem_enum(): assert TwoFlavor[TwoFlavor.NU_E] == TwoFlavor.NU_E @@ -43,7 +48,7 @@ def test_makeFlavorScheme(): TestFlavor = FlavorScheme.from_lepton_names('TestFlavor',leptons=['A','B','C']) assert len(TestFlavor)==6 assert [f.name for f in TestFlavor]==['NU_A','NU_A_BAR','NU_B','NU_B_BAR','NU_C','NU_C_BAR'] - + @staticmethod def test_flavor_properties(): f = ThreeFlavor.NU_E @@ -104,7 +109,12 @@ def test_getitem(): assert m['NU_E','NU_X']==0 assert np.allclose(m['NU_E'], [1,0,0,0]) assert np.allclose(m['NU_E'], m['NU_E',:]) - assert np.allclose(m[:,:], m.array) + + @staticmethod + def test_getitem_submatrix(): + m = FlavorMatrix.eye(TwoFlavor) + assert np.allclose(m[['e','x'],['e','x']], [[1,0],[0,1]]) + assert np.allclose(m[:,:].array, m.array) @staticmethod def test_getitem_short(): From 14a80a6d998e58f5414d0194f01d9556289f0f31 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Wed, 31 Jul 2024 18:02:00 +0300 Subject: [PATCH 41/42] Start working on QuantumDecoherence --- python/snewpy/flavor_transformation.py | 43 +++----------------------- 1 file changed, 5 insertions(+), 38 deletions(-) diff --git a/python/snewpy/flavor_transformation.py b/python/snewpy/flavor_transformation.py index f947d1ee9..435854e4b 100644 --- a/python/snewpy/flavor_transformation.py +++ b/python/snewpy/flavor_transformation.py @@ -418,23 +418,9 @@ def gamma(self, E): return self.m*c.c / (E*self.tau) def get_probabilities(self, t, E): - """neutrino and antineutrino transition probabilities. - - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - PND : an array of 6 x 6 matrices - """ decay_factor = np.exp(-self.gamma(E)*self.d) - PND = FlavorMatrix.eye(self.mix_params.basis_mass, extra_dims=[len(E)]) - + if self.mix_params.mass_order == MassHierarchy.NORMAL: PND['NU_1','NU_3'] = 1 - decay_factor PND['NU_3','NU_3'] = decay_factor @@ -443,7 +429,7 @@ def get_probabilities(self, t, E): PND['NU_2','NU_2'] = decay_factor PND['NU_3','NU_2'] = 1 - decay_factor - PND[3:,3:] = PND[:3,:3] + PND['NU_BAR','NU_BAR'] = PND['NU','NU'] return PND ############################################################################### @@ -482,24 +468,9 @@ def __init__(self, mix_params=None, Gamma3=1e-27*u.eV, Gamma8=1e-27*u.eV, dist=1 self.n = n self.E0 = E0 - def __str__(self): - return f'QuantumDecoherence_' + str(self.mix_params.mass_order) - def get_probabilities(self, t, E): - """neutrino and antineutrino transition probabilities. - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - PQD : an array of 6 x 6 matrices - """ - PQD = np.zeros((6,6,len(E))) + PQD = FlavorMatrix.zeros(self.mix_params.basis_mass, extra_dims=(len(E))) PQD[1,1] = 1/3 + 1/2 * np.exp(-(self.Gamma3 * (E/self.E0)**self.n + self.Gamma8 * (E/self.E0)**self.n / 3) * self.d) \ + 1/6 * np.exp(-self.Gamma8 * (E/self.E0)**self.n * self.d) @@ -509,7 +480,7 @@ def get_probabilities(self, t, E): PQD[1,3] = 1/3 - 1/3 * np.exp(-self.Gamma8 * (E/self.E0)**self.n * self.d) - PQD[2,1] = PQD[1,2] + PQD[2,[1,2,3]] = PQD[1,[2,1,3]] PQD[2,2] = PQD[1,1] PQD[2,3] = PQD[1,3] @@ -517,11 +488,7 @@ def get_probabilities(self, t, E): PQD[3,2] = PQD[2,3] PQD[3,3] = 1/3 + 2/3 * np.exp(-self.Gamma8 * (E/self.E0)**self.n * self.d) - - for i in range(3): - for j in range(3): - PQD[i+3,j+3] = PQD[i,j] - + PQD['NU_BAR','NU_BAR'] = PQD['NU','NU'] return PQD ############################################################################### From 818150348d8258ae0978434c644a0faf25069187 Mon Sep 17 00:00:00 2001 From: Andrey Sheshukov Date: Thu, 1 Aug 2024 11:47:16 +0300 Subject: [PATCH 42/42] Update _repr_markdown_ to prevent creating unnecessary headers --- python/snewpy/flavor.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/python/snewpy/flavor.py b/python/snewpy/flavor.py index 137053ecf..b81e9dc83 100644 --- a/python/snewpy/flavor.py +++ b/python/snewpy/flavor.py @@ -148,13 +148,11 @@ def _repr_markdown_(self): res = [f'**{self.__class__.__name__}**:<`{self.flavor_in.__name__}`->`{self.flavor_out.__name__}`> shape={self.shape}',''] flavors0 = [f.to_tex() for f in self.flavor_in] flavors1 = [f.to_tex() for f in self.flavor_out] - hdr = '|'.join(['*']+flavors1) - res+=[hdr] + res+=['|'.join(['*']+flavors1)] res+=['|'.join(['-:',]*(len(flavors1)+1))] for f0 in self.flavor_in: line = [f0.to_tex()]+[f'{v:.3g}' for v in self[:,f0]] res+=['|'.join(line)] - res+=['---------'] return '\n'.join(res) def __eq__(self,other):