diff --git a/README.md b/README.md index 2aad5c2a..2c43f245 100644 --- a/README.md +++ b/README.md @@ -23,6 +23,17 @@ Alternatively, you can run the following command to explicitly download models y `python -c 'import snewpy; snewpy.get_models()'` +### Earth-Matter Effect +The EMEWS module is available [here](https://github.com/SNEWS2/EMEWS). Follow the installation instructions provided with the code. Whenever you use the module you need to provide a density and an electron fraction profile of the Earth. The PREM is provided. The script SNEWS2.0_rate_table_singleexample+EME.py found in python/snewpy/scripts uses the module. + + +### For Developers + +**Your contributions to SNEWPY are welcome!** For minor changes, simply submit a pull request. If you plan larger changes, it’s probably a good idea to open an issue first to coordinate our work. + +To contribute, first clone the repository (`git clone https://github.com/SNEWS2/snewpy.git`), then make changes and install your modified version locally using `pip install .` from the base directory of the repository. +Once you’re happy with your changes, please submit a pull request. +Unit tests will run automatically for every pull request or you can run them locally using `python -m unittest python/snewpy/test/test_*.py`. ## Usage and Documentation @@ -66,4 +77,4 @@ For more, see the [full documentation on Read the Docs](https://snewpy.rtfd.io/) **Your contributions to SNEWPY are welcome!** For minor changes, simply submit a pull request. If you plan larger changes, it’s probably a good idea to open an issue first to coordinate our work. We use a [Fork & Pull Request](https://docs.github.com/en/get-started/quickstart/fork-a-repo) workflow, which is common on GitHub. -Please see the [Contributing page](https://snewpy.readthedocs.io/en/stable/contributing.html) in our full documentation for details. \ No newline at end of file +Please see the [Contributing page](https://snewpy.readthedocs.io/en/stable/contributing.html) in our full documentation for details. diff --git a/python/snewpy/flavor_transformation.py b/python/snewpy/flavor_transformation.py index d30959ce..18381bb8 100644 --- a/python/snewpy/flavor_transformation.py +++ b/python/snewpy/flavor_transformation.py @@ -10,14 +10,35 @@ 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, MixingParameters +from .neutrino import MassHierarchy, Flavor +try: + import EMEWS +except: + EMEWS = None 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 : 4x4 array or array of 4 x 4 arrays + """ + pass + def prob_ee(self, t, E): """Electron neutrino survival probability. @@ -33,11 +54,11 @@ def prob_ee(self, t, E): float or ndarray Transition probability. """ - pass + p = get_probabilities(t, E) + return p[Flavor.NU_E,Flavor.NU_E] - @abstractmethod def prob_ex(self, t, E): - """X -> e neutrino transition probability. + """x -> e neutrino transition probability. Parameters ---------- @@ -51,11 +72,11 @@ def prob_ex(self, t, E): float or ndarray Transition probability. """ - pass + p = get_probabilities(t, E) + return p[Flavor.NU_E,Flavor.NU_X] - @abstractmethod def prob_xx(self, t, E): - """Flavor X neutrino survival probability. + """Flavor x neutrino survival probability. Parameters ---------- @@ -69,11 +90,11 @@ def prob_xx(self, t, E): float or ndarray Transition probability. """ - pass + p = get_probabilities(t, E) + return p[Flavor.NU_X,Flavor.NU_X] - @abstractmethod def prob_xe(self, t, E): - """e -> X neutrino transition probability. + """e -> x neutrino transition probability. Parameters ---------- @@ -87,9 +108,9 @@ def prob_xe(self, t, E): float or ndarray Transition probability. """ - pass + p = get_probabilities(t, E) + return p[Flavor.NU_X,Flavor.NU_E] - @abstractmethod def prob_eebar(self, t, E): """Electron antineutrino survival probability. @@ -105,11 +126,11 @@ def prob_eebar(self, t, E): float or ndarray Transition probability. """ - pass + p = get_probabilities(t, E) + return p[Flavor.NU_E_BAR,Flavor.NU_E_BAR] - @abstractmethod def prob_exbar(self, t, E): - """X -> e antineutrino transition probability. + """x -> e antineutrino transition probability. Parameters ---------- @@ -123,11 +144,11 @@ def prob_exbar(self, t, E): float or ndarray Transition probability. """ - pass + p = get_probabilities(t, E) + return p[Flavor.NU_E_BAR,Flavor.NU_X_BAR] - @abstractmethod def prob_xxbar(self, t, E): - """X -> X antineutrino survival probability. + """x -> x antineutrino survival probability. Parameters ---------- @@ -141,11 +162,11 @@ def prob_xxbar(self, t, E): float or ndarray Transition probability. """ - pass + p = get_probabilities(t, E) + return p[Flavor.NU_X_BAR,Flavor.NU_X_BAR] - @abstractmethod def prob_xebar(self, t, E): - """e -> X antineutrino transition probability. + """e -> x antineutrino transition probability. Parameters ---------- @@ -159,136 +180,227 @@ def prob_xebar(self, t, E): float or ndarray Transition probability. """ - pass - + p = get_probabilities(t, E) + return p[Flavor.NU_X_BAR,Flavor.NU_E_BAR] -class NoTransformation(FlavorTransformation): - """Survival probabilities for no oscillation case.""" - - def __init__(self): - pass - def prob_ee(self, t, E): - """Electron neutrino survival probability. +class ThreeFlavorTransformation(FlavorTransformation): + """Base class defining common data and methods for all three flavor transformations""" + def __init__(self, mix_params = None, AltAz = None): + """Initialize flavor transformation + Parameters ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - - Returns - ------- - prob : float or ndarray - Transition probability. + mix_params : MixingParameters3Flavor instance or None + AltAz : astropy AltAz object + If not None the Earth-matter effect calculation will be done """ - return 1. + if mix_params == None: + self.mix_params = MixingParameters3Flavor(MassHierarchy.NORMAL) + else: + self.mix_params = mix_params + + self.AltAz = AltAz - def prob_ex(self, t, E): - """X -> e neutrino transition probability. + self.DV_e1 = float((np.cos(mix_params.theta12) * np.cos(mix_params.theta13))**2) + self.DV_e2 = float((np.sin(mix_params.theta12) * np.cos(mix_params.theta13))**2) + self.DV_e3 = float((np.sin(mix_params.theta13))**2) - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. + self.D_e1 = None + self.D_e2 = None + self.D_e3 = None + self.Dbar_e1 = None + self.Dbar_e2 = None + self.Dbar_e3 = None - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_ee(t,E) + self.prior_E = None - def prob_xx(self, t, E): - """Flavor X neutrino survival probability. + + def vacuum(self, E): + """the first rows of the D and Dbar matrices for the case of no Earth matter effect Parameters ---------- - t : float or ndarray - List of times. E : float or ndarray List of energies. Returns ------- - prob : float or ndarray - Transition probability. + the six floats or six ndarrays of the first rows of the D and Dbar matrices """ - return (1. + self.prob_ee(t,E)) / 2. + if self.prior_E != None: + if u.allclose(self.prior_E, E) == True: + return self.D_e1, self.D_e2, self.D_e3, self.Dbar_e1, self.Dbar_e2, self.Dbar_e3 - def prob_xe(self, t, E): - """e -> X neutrino transition probability. + if self.prior_E == None: + self.D_e1 = np.zeros(len(E)) + self.D_e2 = np.zeros(len(E)) + self.D_e3 = np.zeros(len(E)) + self.Dbar_e1 = np.zeros(len(E)) + self.Dbar_e2 = np.zeros(len(E)) + self.Dbar_e3 = np.zeros(len(E)) - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. + self.prior_E = E - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_ee(t,E)) / 2. + self.D_e1 = np.full(len(E),self.DV_e1) + self.D_e2 = np.full(len(E),self.DV_e2) + self.D_e3 = np.full(len(E),self.DV_e3) - def prob_eebar(self, t, E): - """Electron antineutrino survival probability. + self.Dbar_e1 = np.full(len(E),self.DV_e1) + self.Dbar_e2 = np.full(len(E),self.DV_e2) + self.Dbar_e3 = np.full(len(E),self.DV_e3) + + return self.D_e1, self.D_e2, self.D_e3, self.Dbar_e1, self.Dbar_e2, self.Dbar_e3 + + + def Earth_matter(self, 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 ------- - prob : float or ndarray - Transition probability. + an array of length of the E array with each element being the D matrix after computing Earth-matter effects """ - return 1. + if EMEWS == None: + print("The EMEWS module cannot be found. Results do not include the Earth-matter effect.") + return self.vacuum(E) - def prob_exbar(self, t, E): - """X -> e antineutrino transition probability. + if self.prior_E != None: + if u.allclose(self.prior_E, E) == True: + return self.D + + self.D_e1 = np.zeros(len(E)) + self.D_e2 = np.zeros(len(E)) + self.D_e3 = np.zeros(len(E)) + self.Dbar_e1 = np.zeros(len(E)) + self.Dbar_e2 = np.zeros(len(E)) + self.Dbar_e3 = np.zeros(len(E)) + self.prior_E = E + + E = E.to_value('MeV') + + #input data object for EMEWS + ID = EMEWS.InputDataEMEWS() + + ID.altitude = self.AltAz.alt.deg + ID.azimuth = self.AltAz.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 + 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 indici + Pfm = EMEWS.Run(ID) + + m = 0 + while m X antineutrino survival probability. + self.DV_e1 = float((np.cos(mix_params.theta12) * np.cos(mix_params.theta13) * np.cos(mix_params.theta14))**2) + self.DV_e2 = float((np.sin(mix_params.theta12) * np.cos(mix_params.theta13) * np.cos(mix_params.theta14))**2) + self.DV_e3 = float((np.sin(mix_params.theta13) * np.cos(mix_params.theta14))**2) + self.DV_e4 = float((np.sin(mix_params.theta14))**2) + + self.DV_s1 = float((np.cos(mix_params.theta12) * np.cos(mix_params.theta13) * np.sin(mix_params.theta14))**2) + self.DV_s2 = float((np.sin(mix_params.theta12) * np.cos(mix_params.theta13) * np.sin(mix_params.theta14))**2) + self.DV_s3 = float((np.sin(mix_params.theta13) * np.sin(mix_params.theta14))**2) + self.DV_s4 = float((np.cos(mix_params.theta14))**2) + + + def vacuum(self, E): + """the first and last rows of the D and Dbar matrices for the case of no Earth matter effect Parameters ---------- - t : float or ndarray - List of times. E : float or ndarray List of energies. Returns ------- - prob : float or ndarray - Transition probability. + the 16 floats or 16 arrays of the first and last rows of the D and Dbar matrices """ - return (1. + self.prob_eebar(t,E)) / 2. + D_e1 = np.full(len(E),self.DV_e1) + D_e2 = np.full(len(E),self.DV_e2) + D_e3 = np.full(len(E),self.DV_e3) + D_e4 = np.full(len(E),self.DV_e4) - def prob_xebar(self, t, E): - """e -> X antineutrino transition probability. + D_s1 = np.full(len(E),self.DV_s1) + D_s2 = np.full(len(E),self.DV_s2) + D_s3 = np.full(len(E),self.DV_s3) + D_s4 = np.full(len(E),self.DV_s4) + + Dbar_e1 = np.full(len(E),self.DV_e1) + Dbar_e2 = np.full(len(E),self.DV_e2) + Dbar_e3 = np.full(len(E),self.DV_e3) + Dbar_e4 = np.full(len(E),self.DV_e4) + + Dbar_s1 = np.full(len(E),self.DV_s1) + Dbar_s2 = np.full(len(E),self.DV_s2) + Dbar_s3 = np.full(len(E),self.DV_s3) + Dbar_s4 = np.full(len(E),self.DV_s4) + + return D_e1, D_e2, D_e3, D_e4, D_s1, D_s2, D_s3, D_s4, Dbar_e1, Dbar_e2, Dbar_e3, Dbar_e4, Dbar_s1, Dbar_s2, Dbar_s3, Dbar_s4 + + +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 ---------- @@ -299,10 +411,24 @@ def prob_xebar(self, t, E): Returns ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_eebar(t,E)) / 2. + p : 4x4 array or array of 4 x 4 arrays + """ + p_ee = np.ones(len(E)) #p_ee is the probability the an electron neutrino will stay an electron neutrino + pbar_ee = np.ones(len(E)) #pee_bar is the probability that an electron antineutrino will stay an electron antineutrino + + p = np.zeros((4,4,len(E))) + + p[Flavor.NU_E,Flavor.NU_E] = p_ee + p[Flavor.NU_E,Flavor.NU_X] = 1 - p_ee + p[Flavor.NU_X,Flavor.NU_E] = (1 - p_ee)/2 + p[Flavor.NU_X,Flavor.NU_X] = (1 + p_ee)/2 + + p[Flavor.NU_E_BAR,Flavor.NU_E_BAR] = pbar_ee + p[Flavor.NU_E_BAR,Flavor.NU_X_BAR] = 1 - pbar_ee + p[Flavor.NU_X_BAR,Flavor.NU_E_BAR] = (1 - pbar_ee)/2 + p[Flavor.NU_X_BAR,Flavor.NU_X_BAR] = (1 + pbar_ee)/2 + + return p class CompleteExchange(FlavorTransformation): @@ -311,157 +437,171 @@ class CompleteExchange(FlavorTransformation): 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 __str__(self): + return f'CompleteExchange' - 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 get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities. - 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. + p : 4x4 array or array of 4 x 4 arrays + """ + p_ee = np.zeros(len(E)) #p_ee is the probability the an electron neutrino will stay an electron neutrino + pbar_ee = np.zeros(len(E)) #pee_bar is the probability that an electron antineutrino will stay an electron antineutrino - def prob_eebar(self, t, E): - """Electron antineutrino survival probability. + p = np.zeros((4,4,len(E))) + + p[Flavor.NU_E,Flavor.NU_E] = p_ee + p[Flavor.NU_E,Flavor.NU_X] = 1 - p_ee + p[Flavor.NU_X,Flavor.NU_E] = (1 - p_ee)/2 + p[Flavor.NU_X,Flavor.NU_X] = (1 + p_ee)/2 + + p[Flavor.NU_E_BAR,Flavor.NU_E_BAR] = pbar_ee + p[Flavor.NU_E_BAR,Flavor.NU_X_BAR] = 1 - pbar_ee + p[Flavor.NU_X_BAR,Flavor.NU_E_BAR] = (1 - pbar_ee)/2 + p[Flavor.NU_X_BAR,Flavor.NU_X_BAR] = (1 + pbar_ee)/2 + + return p + + +class AdiabaticMSW(ThreeFlavorTransformation): + """Adiabatic MSW effect.""" + + def __init__(self, mix_params, AltAz=None): + """Initialize flavor transformation + Parameters ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - Returns - ------- - prob : float or ndarray - Transition probability. + mix_params : MixingParameters3Flavor instance + AltAz : astropy AltAz object + If not None the Earth-matter effect calculation will be done """ - return 0. + super().__init__(mix_params, AltAz) + + def __str__(self): + return f'AdiabaticMSW' + + def get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities. - 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) + p : 4x4 array or array of 4 x 4 arrays + """ + if self.AltAz == None or self.AltAz.alt > 0: + D_e1, D_e2, D_e3, Dbar_e1, Dbar_e2, Dbar_e3 = self.vacuum(E) + else: + D_e1, D_e2, D_e3, Dbar_e1, Dbar_e2, Dbar_e3 = self.Earth_matter(E) - def prob_xxbar(self, t, E): - """X -> X antineutrino survival probability. + if self.mix_params.mass_order == MassHierarchy.NORMAL: + p_ee = D_e3 + pbar_ee = Dbar_e1 + else: + p_ee = D_e2 + pbar_ee = Dbar_e3 + + p = np.zeros((4,4,len(E))) + + p[Flavor.NU_E,Flavor.NU_E] = p_ee + p[Flavor.NU_E,Flavor.NU_X] = 1 - p_ee + p[Flavor.NU_X,Flavor.NU_E] = (1 - p_ee)/2 + p[Flavor.NU_X,Flavor.NU_X] = (1 + p_ee)/2 + + p[Flavor.NU_E_BAR,Flavor.NU_E_BAR] = pbar_ee + p[Flavor.NU_E_BAR,Flavor.NU_X_BAR] = 1 - pbar_ee + p[Flavor.NU_X_BAR,Flavor.NU_E_BAR] = (1 - pbar_ee)/2 + p[Flavor.NU_X_BAR,Flavor.NU_X_BAR] = (1 + pbar_ee)/2 + + return p + + +class NonAdiabaticMSWH(ThreeFlavorTransformation): + """Nonadiabatic MSW effect.""" + + def __init__(self, mix_params, AltAz=None): + """Initialize flavor transformation + Parameters ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. - Returns - ------- - prob : float or ndarray - Transition probability. + mix_params : MixingParameters3Flavor instance + AltAz : astropy AltAz object + If not None the Earth-matter effect calculation will be done """ - return (1. + self.prob_eebar(t,E)) / 2. + super().__init__(mix_params, AltAz) + + def __str__(self): + return f'NonAdiabaticMSWH' + + def get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities. - 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. + if self.AltAz == None or self.AltAz.alt > 0: + D_e1, D_e2, D_e3, Dbar_e1, Dbar_e2, Dbar_e3 = self.vacuum(E) + else: + D_e1, D_e2, D_e3, Dbar_e1, Dbar_e2, Dbar_e3 = self.Earth_matter(E) + + if self.mix_params.mass_order == MassHierarchy.NORMAL: + p_ee = D_e2 + pbar_ee = Dbar_e1 + else: + p_ee = D_e2 + pbar_ee = Dbar_e1 + + p = np.zeros((4,4,len(E))) + T[Flavor.NU_E,Flavor.NU_E] = p_ee + p[Flavor.NU_E,Flavor.NU_X] = 1 - p_ee + p[Flavor.NU_X,Flavor.NU_E] = (1 - p_ee)/2 + p[Flavor.NU_X,Flavor.NU_X] = (1 + p_ee)/2 -class AdiabaticMSW(FlavorTransformation): - """Adiabatic MSW effect.""" + p[Flavor.NU_E_BAR,Flavor.NU_E_BAR] = pbar_ee + p[Flavor.NU_E_BAR,Flavor.NU_X_BAR] = 1 - pbar_ee + p[Flavor.NU_X_BAR,Flavor.NU_E_BAR] = (1 - pbar_ee)/2 + p[Flavor.NU_X_BAR,Flavor.NU_X_BAR] = (1 + pbar_ee)/2 - def __init__(self, mix_angles=None, mh=MassHierarchy.NORMAL): - """Initialize transformation matrix. + return p + +class TwoFlavorDecoherence(ThreeFlavorTransformation): + """Star-earth transit survival probability: two flavor case.""" + def __init__(self, mix_params, AltAz=None): + """Initialize flavor transformation + 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. + mix_params : MixingParameters3Flavor instance + AltAz : astropy AltAz object + If not None the Earth-matter effect calculation will be done """ - 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() + super().__init__(mix_params, AltAz) - 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 __str__(self): + return f'TwoFlavorDecoherence' - def prob_ee(self, t, E): - """Electron neutrino survival probability. + def get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities. Parameters ---------- @@ -472,33 +612,46 @@ def prob_ee(self, t, E): Returns ------- - prob : float or ndarray - Transition probability. - """ - if self.mass_order == MassHierarchy.NORMAL: - return self.De3 + p : 4x4 array or array of 4 x 4 arrays + """ + if self.AltAz == None or self.AltAz.alt > 0: + D_e1, D_e2, D_e3, Dbar_e1, Dbar_e2, Dbar_e3 = self.vacuum(E) else: - return self.De2 + D_e1, D_e2, D_e3, Dbar_e1, Dbar_e2, Dbar_e3 = self.Earth_matter(E) - def prob_ex(self, t, E): - """X -> e neutrino transition probability. + if self.mix_params.mass_order == MassHierarchy.NORMAL: + p_ee = (D_e2 + D_e3)/2 + pbar_ee = Dbar_e1 + else: + p_ee = D_e2 + pbar_ee = (Dbar_e1 + Dbar_e3)/2 - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. + p = np.zeros((4,4,len(E))) - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_ee(t,E) + p[Flavor.NU_E,Flavor.NU_E] = p_ee + p[Flavor.NU_E,Flavor.NU_X] = 1 - p_ee + p[Flavor.NU_X,Flavor.NU_E] = (1 - p_ee)/2 + p[Flavor.NU_X,Flavor.NU_X] = (1 + p_ee)/2 - def prob_xx(self, t, E): - """Flavor X neutrino survival probability. + p[Flavor.NU_E_BAR,Flavor.NU_E_BAR] = pbar_ee + p[Flavor.NU_E_BAR,Flavor.NU_X_BAR] = 1 - pbar_ee + p[Flavor.NU_X_BAR,Flavor.NU_E_BAR] = (1 - pbar_ee)/2 + p[Flavor.NU_X_BAR,Flavor.NU_X_BAR] = (1 + pbar_ee)/2 + + return p + + +class ThreeFlavorDecoherence(FlavorTransformation): + """Star-earth transit survival probability: three flavor case.""" + + def __init__(self): + pass + + def __str__(self): + return f'ThreeFlavorDecoherence' + + def get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities. Parameters ---------- @@ -509,751 +662,76 @@ def prob_xx(self, t, E): Returns ------- - prob : float or ndarray - Transition probability. - """ - return (1. + self.prob_ee(t,E)) / 2. + p : 4x4 array or array of 4 x 4 arrays + """ + p_ee.full(len(E),1/3) + pbar_ee.full(len(E),1/3) - def prob_xe(self, t, E): - """e -> X neutrino transition probability. + p = np.zeros((4,4,len(E))) - 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. + p[Flavor.NU_E,Flavor.NU_E] = p_ee + p[Flavor.NU_E,Flavor.NU_X] = 1 - p_ee + p[Flavor.NU_X,Flavor.NU_E] = (1 - p_ee)/2 + p[Flavor.NU_X,Flavor.NU_X] = (1 + p_ee)/2 - 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. + p[Flavor.NU_E_BAR,Flavor.NU_E_BAR] = pbar_ee + p[Flavor.NU_E_BAR,Flavor.NU_X_BAR] = 1 - pbar_ee + p[Flavor.NU_X_BAR,Flavor.NU_E_BAR] = (1 - pbar_ee)/2 + p[Flavor.NU_X_BAR,Flavor.NU_X_BAR] = (1 + pbar_ee)/2 - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_eebar(t,E)) / 2. + return p -class NeutrinoDecay(FlavorTransformation): +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_angles=None, mass=1*u.eV/c.c**2, tau=1*u.day, dist=10*u.kpc, mh=MassHierarchy.NORMAL): + def __init__(self, mix_params, AltAz=None, mass=1*u.eV/c.c**2, tau=1*u.day, dist=10*u.kpc): """Initialize transformation matrix. Parameters ---------- - mix_angles : tuple or None - If not None, override default mixing angles using tuple (theta12, theta13, theta23). + mix_params : MixingParameters3Flavor instance + AltAz : astropy AltAz object + If not None the Earth-matter effect calculation will be done 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. + 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. + AltAz : astropy AltAz object + If not None the Earth-matter effect calculation will be done """ - return self.De3 + super().__init__(mix_params, AltAz) + + self.m = mass + self.tau = tau + self.d = dist - def prob_exbar(self, t, E): - """X -> e antineutrino transition probability. + def __str__(self): + return f'NeutrinoDecay' + + def gamma(self, E): + """Decay width of the heaviest neutrino mass. Parameters ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. + E : float + Energy of the nu3. Returns ------- - prob : float or ndarray - Transition probability. + Gamma : float + Decay width of the neutrino mass, in units of 1/length. + + :meta private: """ - # 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 + return self.m*c.c / (E*self.tau) - def prob_xxbar(self, t, E): - """X -> X antineutrino survival probability. + def get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities. Parameters ---------- @@ -1264,30 +742,42 @@ def prob_xxbar(self, t, E): Returns ------- - prob : float or ndarray - Transition probability. - """ - return 1. - self.prob_exbar(t,E) / 2. + p : 4x4 array or array of 4 x 4 arrays + """ + if self.AltAz == None or self.AltAz.alt > 0: + D_e1, D_e2, D_e3, Dbar_e1, Dbar_e2, Dbar_e3 = self.vacuum(E) + else: + D_e1, D_e2, D_e3, Dbar_e1, Dbar_e2, Dbar_e3 = self.Earth_matter(E) + + decay_factor = math.e**(-self.gamma(E)*self.dist) + + if self.mix_params.mass_order == MassHierarchy.NORMAL: + p_ee = D_e1 * (1 - decay_factor) + D_e3 * decay_factor + p_ex = D_e1 + D_e2 + pbar_ee = Dbar_e1 + pbar_ex = Dbar_e1 * (1 - decay_factor) + Dbar_e2 + Dbar_e3 * decay_factor + else: + p_ee = D_e2 * decay_factor + D_e3 * (1 - decay_factor) + p_ex = D_e1 + D_e3 + pbar_ee = Dbar_e3 + pbar_ex = Dbar_e1 + Dbar_e2 * decay_factor + Dbar_e3 * (1 - decay_factor) - def prob_xebar(self, t, E): - """e -> X antineutrino transition probability. + p = np.zeros((4,4,len(E))) - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. + p[Flavor.NU_E,Flavor.NU_E] = p_ee + p[Flavor.NU_E,Flavor.NU_X] = p_ex + p[Flavor.NU_X,Flavor.NU_E] = (1 - p_ee)/2 + p[Flavor.NU_X,Flavor.NU_X] = 1 - p_ex/2 - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return (1. - self.prob_eebar(t,E)) / 2. + p[Flavor.NU_E_BAR,Flavor.NU_E_BAR] = pbar_ee + p[Flavor.NU_E_BAR,Flavor.NU_X_BAR] = pbar_ex + p[Flavor.NU_X_BAR,Flavor.NU_E_BAR] = (1 - pbar_ee)/2 + p[Flavor.NU_X_BAR,Flavor.NU_X_BAR] = 1 - pbar_ex/2 + + return p -class AdiabaticMSWes(FlavorTransformation): +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 @@ -1297,108 +787,20 @@ class AdiabaticMSWes(FlavorTransformation): 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 __init__(self, mix_params): + """Initialize flavor transformation - 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. + mix_params : MixingParameters3Flavor instance + """ + super().__init__(mix_params) - Returns - ------- - prob : float or ndarray - Transition probability. - """ - return ( 1 - self.De4 - self.Ds4 )/2 + def __str__(self): + return f'AdiabaticMSWes' - def prob_eebar(self, t, E): - """e -> e antineutrino transition probability. + def get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities. Parameters ---------- @@ -1409,76 +811,47 @@ def prob_eebar(self, t, E): 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. + p : 4x4 array or array of 4 x 4 arrays + """ + D_e1, D_e2, D_e3, D_e4, D_s1, D_s2, D_s3, D_s4, Dbar_e1, Dbar_e2, Dbar_e3, Dbar_e4, Dbar_s1, Dbar_s2, Dbar_s3, Dbar_s4 = self.vacuum(E) - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. + if self.mix_params.mass_order == MassHierarchy.NORMAL: + p_ee = D_e4 + p_ex = D_e1 + D_e2 + p_xe = (1 - D_e4 - D_s4) / 2 + p_xx = (2 - D_e1 - D_e2 - D_s1 - D_s2) / 2 - Returns - ------- - prob : float or ndarray - Transition probability. - """ - if self.mass_order == MassHierarchy.NORMAL: - return self.De3 + self.De4 + pbar_ee = Dbar_e1 + pbar_ex = Dbar_e3 + Dbar_e4 + pbar_xe = (1 - Dbar_e1 - Dbar_s1) / 2 + pbar_xx = (2 - Dbar_e3 - Dbar_e4 - Dbar_s3 - Dbar_s4) / 2 else: - return self.De2 + self.De4 - - def prob_xxbar(self, t, E): - """x -> x antineutrino transition probability. + p_ee = D_e4 + p_ex = D_e1 + D_e3 + p_xe = (1 - D_e4 - D_s4) / 2 + p_xx = (2 - D_e1 - D_e3 - D_s1 - D_s3) / 2 - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. + pbar_ee = Dbar_e4 + pbar_ex = Dbar_e1 + Dbar_e3 + pbar_xe = (1 - Dbar_e4 - Dbar_s4) / 2 + pbar_xx = (2 - Dbar_e1 - Dbar_e3 - Dbar_s1 - Dbar_s3) / 2 - 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. + p = np.zeros((4,4,len(E))) - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. + p[Flavor.NU_E,Flavor.NU_E] = p_ee + p[Flavor.NU_E,Flavor.NU_X] = p_ex + p[Flavor.NU_X,Flavor.NU_E] = p_xe + p[Flavor.NU_X,Flavor.NU_X] = p_xx - 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 + p[Flavor.NU_E_BAR,Flavor.NU_E_BAR] = pbar_ee + p[Flavor.NU_E_BAR,Flavor.NU_X_BAR] = pbar_ex + p[Flavor.NU_X_BAR,Flavor.NU_E_BAR] = pbar_xe + p[Flavor.NU_X_BAR,Flavor.NU_X_BAR] = pbar_xx + return p + -class NonAdiabaticMSWes(FlavorTransformation): +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 @@ -1488,94 +861,20 @@ class NonAdiabaticMSWes(FlavorTransformation): 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. + def __init__(self, mix_params, mh=MassHierarchy.NORMAL): + """Initialize flavor transformation Parameters ---------- - mix_angles : tuple - Values for mixing angles (theta12, theta13, theta23, theta14). - mh : MassHierarchy - MassHierarchy.NORMAL or MassHierarchy.INVERTED. + mix_params : MixingParameters3Flavor instance """ - 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. + super().__init__(mix_params) - 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. + def __str__(self): + return f'NonAdiabaticMSWes' + + def get_probabilities(self, t, E): + """neutrino and antineutrino transition probabilities. Parameters ---------- @@ -1586,91 +885,43 @@ def prob_xe(self, t, E): Returns ------- - prob : float or ndarray - Transition probability. + p : 4x4 array or array of 4 x 4 arrays """ - if self.mass_order == MassHierarchy.NORMAL: - return (1 - self.De3 - self.Ds3)/2 + D_e1, D_e2, D_e3, D_e4, D_s1, D_s2, D_s3, D_s4, Dbar_e1, Dbar_e2, Dbar_e3, Dbar_e4, Dbar_s1, Dbar_s2, Dbar_s3, Dbar_s4 = self.vacuum(E) + + if self.mix_params.mass_order == MassHierarchy.NORMAL: + p_ee = D_e3 + p_ex = D_e1 + D_e2 + p_xe = (1 - D_e3) / 2 + p_xx = (2 - D_e1 - D_e2 - D_s1 - D_s2) / 2 + + pbar_ee = Dbar_e1 + pbar_ex = Dbar_e2 + Dbar_e3 + pbar_xe = (1 - Dbar_e1) / 2 + pbar_xx = (2 - Dbar_e2 - Dbar_e3 - Dbar_s2 - Dbar_s3) / 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. + p_ee = D_e2 + p_ex = D_e1 + D_e3 + p_xe = (1 - D_e2) / 2 + p_xx = (2 - D_e1 - D_e3 - D_s1 - D_s3) / 2 - 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. + pbar_ee = Dbar_e3 + pbar_ex = Dbar_e1 + Dbar_e2 + pbar_xe = (1 - Dbar_e3) / 2 + pbar_xx = (2 - Dbar_e1 - Dbar_e2 - Dbar_s1 - Dbar_s2) / 2 - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. + p = np.zeros((4,4,len(E))) - 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. + p[Flavor.NU_E,Flavor.NU_E] = p_ee + p[Flavor.NU_E,Flavor.NU_X] = p_ex + p[Flavor.NU_X,Flavor.NU_E] = p_xe + p[Flavor.NU_X,Flavor.NU_X] = p_xx - Parameters - ---------- - t : float or ndarray - List of times. - E : float or ndarray - List of energies. + p[Flavor.NU_E_BAR,Flavor.NU_E_BAR] = pbar_ee + p[Flavor.NU_E_BAR,Flavor.NU_X_BAR] = pbar_ex + p[Flavor.NU_X_BAR,Flavor.NU_E_BAR] = pbar_xe + p[Flavor.NU_X_BAR,Flavor.NU_X_BAR] = pbar_xx - 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 + return p - 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 - + diff --git a/python/snewpy/models/base.py b/python/snewpy/models/base.py index 97f71f76..e2ea7ad1 100644 --- a/python/snewpy/models/base.py +++ b/python/snewpy/models/base.py @@ -145,23 +145,24 @@ def get_transformed_spectra(self, t, E, flavor_xform): Dictionary of transformed spectra, keyed by neutrino flavor. """ initialspectra = self.get_initial_spectra(t, E) + transformation_matrix = flavor_xform.get_probabilities(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] + transformation_matrix[Flavor.NU_E, Flavor.NU_E] * initialspectra[Flavor.NU_E] + \ + transformation_matrix[Flavor.NU_E, Flavor.NU_X] * 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] + transformation_matrix[Flavor.NU_X, Flavor.NU_E] * initialspectra[Flavor.NU_E] + \ + transformation_matrix[Flavor.NU_X, Flavor.NU_X] * 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] - + transformation_matrix[Flavor.NU_E_BAR, Flavor.NU_E_BAR] * initialspectra[Flavor.NU_E_BAR] + \ + transformation_matrix[Flavor.NU_E_BAR, Flavor.NU_X_BAR] * 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] + transformation_matrix[Flavor.NU_X_BAR, Flavor.NU_E_BAR] * initialspectra[Flavor.NU_E_BAR] + \ + transformation_matrix[Flavor.NU_X_BAR, Flavor.NU_X_BAR] * initialspectra[Flavor.NU_X_BAR] return transformed_spectra diff --git a/python/snewpy/scripts/PREM.Ye.dat b/python/snewpy/scripts/PREM.Ye.dat new file mode 100644 index 00000000..808c79c0 --- /dev/null +++ b/python/snewpy/scripts/PREM.Ye.dat @@ -0,0 +1,49 @@ +0.000000E+00 0.466 +2.000000E+07 0.466 +4.000000E+07 0.466 +6.000000E+07 0.466 +8.000000E+07 0.466 +1.000000E+08 0.466 +1.200000E+08 0.466 +1.221500E+08 0.466 +1.221500E+08 0.466 +1.400000E+08 0.466 +1.600000E+08 0.466 +1.800000E+08 0.466 +2.000000E+08 0.466 +2.200000E+08 0.466 +2.400000E+08 0.466 +2.600000E+08 0.466 +2.800000E+08 0.466 +3.000000E+08 0.466 +3.200000E+08 0.466 +3.400000E+08 0.466 +3.480000E+08 0.5 +3.480000E+08 0.5 +3.600000E+08 0.5 +3.630000E+08 0.5 +3.800000E+08 0.5 +4.000000E+08 0.5 +4.200000E+08 0.5 +4.400000E+08 0.5 +4.600000E+08 0.5 +4.800000E+08 0.5 +5.000000E+08 0.5 +5.200000E+08 0.5 +5.400000E+08 0.5 +5.600000E+08 0.5 +5.701000E+08 0.5 +5.701000E+08 0.5 +5.771000E+08 0.5 +5.871000E+08 0.5 +5.971000E+08 0.5 +5.971000E+08 0.5 +6.061000E+08 0.5 +6.151000E+08 0.5 +6.221000E+08 0.5 +6.291000E+08 0.5 +6.346600E+08 0.5 +6.346600E+08 0.5 +6.356000E+08 0.5 +6.356000E+08 0.5 +6.368000E+08 0.5 diff --git a/python/snewpy/scripts/PREM.rho.dat b/python/snewpy/scripts/PREM.rho.dat new file mode 100644 index 00000000..197103c0 --- /dev/null +++ b/python/snewpy/scripts/PREM.rho.dat @@ -0,0 +1,49 @@ +0.000000E+00 1.308848E+01 +2.000000E+07 1.307977E+01 +4.000000E+07 1.305364E+01 +6.000000E+07 1.301009E+01 +8.000000E+07 1.294912E+01 +1.000000E+08 1.287073E+01 +1.200000E+08 1.277493E+01 +1.221500E+08 1.276360E+01 +1.221500E+08 1.216634E+01 +1.400000E+08 1.206924E+01 +1.600000E+08 1.194682E+01 +1.800000E+08 1.180900E+01 +2.000000E+08 1.165478E+01 +2.200000E+08 1.148311E+01 +2.400000E+08 1.129298E+01 +2.600000E+08 1.108335E+01 +2.800000E+08 1.085321E+01 +3.000000E+08 1.060152E+01 +3.200000E+08 1.032726E+01 +3.400000E+08 1.002940E+01 +3.480000E+08 9.903490E+00 +3.480000E+08 5.566450E+00 +3.600000E+08 5.506420E+00 +3.630000E+08 5.491450E+00 +3.800000E+08 5.406810E+00 +4.000000E+08 5.307240E+00 +4.200000E+08 5.207130E+00 +4.400000E+08 5.105900E+00 +4.600000E+08 5.002990E+00 +4.800000E+08 4.897830E+00 +5.000000E+08 4.789830E+00 +5.200000E+08 4.678440E+00 +5.400000E+08 4.563070E+00 +5.600000E+08 4.443170E+00 +5.701000E+08 4.380710E+00 +5.701000E+08 3.992140E+00 +5.771000E+08 3.975840E+00 +5.871000E+08 3.849800E+00 +5.971000E+08 3.723780E+00 +5.971000E+08 3.543250E+00 +6.061000E+08 3.489510E+00 +6.151000E+08 3.359500E+00 +6.221000E+08 3.367100E+00 +6.291000E+08 3.374710E+00 +6.346600E+08 3.380760E+00 +6.346600E+08 2.900000E+00 +6.356000E+08 2.900000E+00 +6.356000E+08 2.600000E+00 +6.368000E+08 2.600000E+00 diff --git a/python/snewpy/scripts/SNEWS2.0_rate_table_singleexample+EME.py b/python/snewpy/scripts/SNEWS2.0_rate_table_singleexample+EME.py new file mode 100644 index 00000000..413c6684 --- /dev/null +++ b/python/snewpy/scripts/SNEWS2.0_rate_table_singleexample+EME.py @@ -0,0 +1,64 @@ +#!/usr/bin/env python +import numpy as np +from astropy import units as u +from astropy.time import Time +from astropy.coordinates import SkyCoord, EarthLocation, AltAz + +from snewpy import snowglobes +from snewpy.flavor_transformation import * +from snewpy.neutrino import * + +SNOwGLoBES_path = None # change to SNOwGLoBES directory if using a custom detector configuration +SNEWPY_model_dir = "/path/to/snewpy/models/" # directory containing model input files + +# skycoordinates of neutrino source +Betelgeuse = SkyCoord.from_name('Betelgeuse') + +# neutrino detector +SuperK = EarthLocation.of_site('SuperK') + +# time when the supernova occured +# the first time option means the neutrinos traveled through the Earth, the second means they did not +time = Time('2021-5-26 14:14:00') +#time = Time('2021-5-26 14:14:00') - 12*u.hour + +# altaz of supernovae at detector +SNaltaz = Betelgeuse.transform_to(AltAz(obstime=time,location=SuperK)) + +distance = 10 # Supernova distance in kpc +detector = "wc100kt30prct" #SNOwGLoBES detector for water Cerenkov +modeltype = 'Bollig_2016' # Model type from snewpy.models +model = 's11.2c' # Name of model + +mix_params = MixingParameters(MassHierarchy.NORMAL) + +transformation = AdiabaticMSW(mix_params,SNaltaz) # Desired flavor transformation + +# Construct file system path of model file and name of output file +model_path = SNEWPY_model_dir + "/" + modeltype + "/" + model +outfile = modeltype + "_" + model + "_" + str(transformation) + +# Now, do the main work: +print("Generating fluence files ...") +tarredfile = snowglobes.generate_fluence(model_path, modeltype, transformation, distance, outfile) + +print("Simulating detector effects with SNOwGLoBES ...") +snowglobes.simulate(SNOwGLoBES_path, tarredfile, detector_input=detector) + +print("Collating results ...") +tables = snowglobes.collate(SNOwGLoBES_path, tarredfile, skip_plots=True) + + +# Use results to print the number of events in different interaction channels +key = f"Collated_{outfile}_{detector}_events_smeared_weighted.dat" +total_events = 0 +for i, channel in enumerate(tables[key]['header'].split()): + if i == 0: + continue + n_events = sum(tables[key]['data'][i]) + total_events += n_events + print(f"{channel:10}: {n_events:.3f} events") + +#Super-K has 32kT inner volume +print("Total events in Super-K-like detector:",0.32*total_events) + diff --git a/python/snewpy/snowglobes.py b/python/snewpy/snowglobes.py index 53e78cde..be2e5455 100644 --- a/python/snewpy/snowglobes.py +++ b/python/snewpy/snowglobes.py @@ -34,12 +34,43 @@ import snewpy.models from snewpy.flavor_transformation import * -from snewpy.neutrino import MassHierarchy +from snewpy.neutrino import MassHierarchy, MixingParameters from snewpy.rate_calculator import RateCalculator, center from snewpy.flux import Container + logger = logging.getLogger(__name__) -def generate_time_series(model_path, model_type, transformation_type, d, output_filename=None, ntbins=30, deltat=None, snmodel_dict={}): +def get_transformation(flavor_transformation_string): + """A function for idetifying the flavor transformation class from a string + + Paramters + --------- + flavor_transformation_string : string + + Returns + ------- + flavor_transformation_class : a flavor_transformation class + """ + + print(flavor_transformation_string) + + IMO_mix_params = MixingParameters(MassHierarchy.INVERTED) + NMO_mix_params = MixingParameters(MassHierarchy.NORMAL) + + if flavor_transformation_string.find("NeutrinoDecay") is True: + print("Using the default values for Neutrino Decay transformation") + + + # Choose flavor transformation. Use dict to associate the transformation name with its class. + # The default mixing paramaters are the normal hierarchy values + flavor_transformation_dict = {'NoTransformation': NoTransformation(), 'CompleteExchange': CompleteExchange(), 'AdiabaticMSW_NMO': AdiabaticMSW(NMO_mix_params), 'AdiabaticMSW_IMO': AdiabaticMSW(IMO_mix_params), 'NonAdiabaticMSWH_NMO': NonAdiabaticMSWH(NMO_mix_params), 'NonAdiabaticMSWH_IMO': NonAdiabaticMSWH(IMO_mix_params), 'TwoFlavorDecoherence': TwoFlavorDecoherence(NMO_mix_params), 'TwoFlavorDecoherence_NMO': TwoFlavorDecoherence(NMO_mix_params), 'TwoFlavorDecoherence_IMO': TwoFlavorDecoherence(IMO_mix_params), 'ThreeFlavorDecoherence': ThreeFlavorDecoherence(), 'NeutrinoDecay_NMO': NeutrinoDecay(NMO_mix_params), 'NeutrinoDecay_IMO': NeutrinoDecay(IMO_mix_params)} + + flavor_transformation_class = flavor_transformation_dict[flavor_transformation_string] + + return flavor_transformation_class + + +def generate_time_series(model_path, model_type, flavor_transformation, d, output_filename=None, ntbins=30, deltat=None, snmodel_dict={}): """Generate time series files in SNOwGLoBES format. This version will subsample the times in a supernova model, produce energy @@ -51,8 +82,9 @@ def generate_time_series(model_path, model_type, transformation_type, d, output_ Input file containing neutrino flux information from supernova model. model_type : str Format of input file. Matches the name of the corresponding class in :py:mod:`snewpy.models`. - transformation_type : str - Name of flavor transformation. See snewpy.flavor_transformation documentation for possible values. + flavor_transformation : str or flavor transformation class + If a string, the class if found using the get_transformation function. + See snewpy.flavor_transformation documentation for possible values of the string. d : int or float Distance to supernova in kpc. output_filename : str or None @@ -72,8 +104,8 @@ def generate_time_series(model_path, model_type, transformation_type, d, output_ model_class = getattr(snewpy.models.ccsn, model_type) # Choose flavor transformation. Use dict to associate the transformation name with its class. - flavor_transformation_dict = {'NoTransformation': NoTransformation(), 'AdiabaticMSW_NMO': AdiabaticMSW(mh=MassHierarchy.NORMAL), 'AdiabaticMSW_IMO': AdiabaticMSW(mh=MassHierarchy.INVERTED), 'NonAdiabaticMSWH_NMO': NonAdiabaticMSWH(mh=MassHierarchy.NORMAL), 'NonAdiabaticMSWH_IMO': NonAdiabaticMSWH(mh=MassHierarchy.INVERTED), 'TwoFlavorDecoherence': TwoFlavorDecoherence(), 'ThreeFlavorDecoherence': ThreeFlavorDecoherence(), 'NeutrinoDecay_NMO': NeutrinoDecay(mh=MassHierarchy.NORMAL), 'NeutrinoDecay_IMO': NeutrinoDecay(mh=MassHierarchy.INVERTED)} - flavor_transformation = flavor_transformation_dict[transformation_type] + if isinstance(flavor_transformation,str) == True: + flavor_transformation = get_transformation(flavor_transformation) model_dir, model_file = os.path.split(os.path.abspath(model_path)) snmodel = model_class(model_path, **snmodel_dict) @@ -96,11 +128,11 @@ def generate_time_series(model_path, model_type, transformation_type, d, output_ tfname = output_filename + '.npz' else: model_file_root, _ = os.path.splitext(model_file) # strip extension (if present) - tfname = f'{model_file_root}.{transformation_type}.{tmin:.3f},{tmax:.3f},{ntbins:d}-{d:.1f}.npz' + tfname = f'{model_file_root}'+str(flavor_transformation)+f'.{tmin:.3f},{tmax:.3f},{ntbins:d}-{d:.1f}.npz' fluence.save(tfname) return tfname -def generate_fluence(model_path, model_type, transformation_type, d, output_filename=None, tstart=None, tend=None, snmodel_dict={}): +def generate_fluence(model_path, model_type, flavor_transformation, d, output_filename=None, tstart=None, tend=None, snmodel_dict={}): """Generate fluence files in SNOwGLoBES format. This version will subsample the times in a supernova model, produce energy @@ -112,8 +144,9 @@ def generate_fluence(model_path, model_type, transformation_type, d, output_file Input file containing neutrino flux information from supernova model. model_type : str Format of input file. Matches the name of the corresponding class in :py:mod:`snewpy.models`. - transformation_type : str - Name of flavor transformation. See snewpy.flavor_transformation documentation for possible values. + flavor_transformation : str or flavor transformation class + If a string, the class if found using the get_transformation function. + See snewpy.flavor_transformation documentation for possible values of the string. d : int or float Distance to supernova in kpc. output_filename : str or None @@ -133,8 +166,8 @@ def generate_fluence(model_path, model_type, transformation_type, d, output_file model_class = getattr(snewpy.models.ccsn, model_type) # Choose flavor transformation. Use dict to associate the transformation name with its class. - flavor_transformation_dict = {'NoTransformation': NoTransformation(), 'AdiabaticMSW_NMO': AdiabaticMSW(mh=MassHierarchy.NORMAL), 'AdiabaticMSW_IMO': AdiabaticMSW(mh=MassHierarchy.INVERTED), 'NonAdiabaticMSWH_NMO': NonAdiabaticMSWH(mh=MassHierarchy.NORMAL), 'NonAdiabaticMSWH_IMO': NonAdiabaticMSWH(mh=MassHierarchy.INVERTED), 'TwoFlavorDecoherence': TwoFlavorDecoherence(), 'ThreeFlavorDecoherence': ThreeFlavorDecoherence(), 'NeutrinoDecay_NMO': NeutrinoDecay(mh=MassHierarchy.NORMAL), 'NeutrinoDecay_IMO': NeutrinoDecay(mh=MassHierarchy.INVERTED)} - flavor_transformation = flavor_transformation_dict[transformation_type] + if isinstance(flavor_transformation,str) == True: + flavor_transformation = get_transformation(flavor_transformation) model_dir, model_file = os.path.split(os.path.abspath(model_path)) snmodel = model_class(model_path, **snmodel_dict) @@ -165,9 +198,10 @@ def generate_fluence(model_path, model_type, transformation_type, d, output_file tfname = output_filename+'.npz' else: model_file_root, _ = os.path.splitext(model_file) # strip extension (if present) - tfname = f'{model_file_root}.{transformation_type}.{times[0]:.3f},{times[1]:.3f},{len(times)-1:d}-{d:.1f}.npz' + tfname = f'{model_file_root}.'+str(flavor_transformation)+f'.{times[0]:.3f},{times[1]:.3f},{len(times)-1:d}-{d:.1f}.npz' fluence.save(tfname) + return tfname def simulate(SNOwGLoBESdir, tarball_path, detector_input="all", verbose=False, *, detector_effects=True): diff --git a/python/snewpy/test/test_04_xforms.py b/python/snewpy/test/test_04_xforms.py index 676c0504..acdc76f7 100644 --- a/python/snewpy/test/test_04_xforms.py +++ b/python/snewpy/test/test_04_xforms.py @@ -1,7 +1,7 @@ # -*- coding: utf-8 -*- import unittest -from snewpy.flavor_transformation import MassHierarchy, MixingParameters +from snewpy.neutrino import MassHierarchy, MixingParameters from snewpy.flavor_transformation import \ NoTransformation, CompleteExchange, \ AdiabaticMSW, NonAdiabaticMSWH, \