From 9ceedc2fe225835a0fb6b1bf7da9d4bdde9c67de Mon Sep 17 00:00:00 2001 From: mvsantos001 Date: Fri, 14 Jul 2023 08:48:09 -0300 Subject: [PATCH 01/11] Taking NeutrinoDecay class as an example --- python/snewpy/flavor_transformation.py | 215 +++++++++++++++++++++++++ 1 file changed, 215 insertions(+) diff --git a/python/snewpy/flavor_transformation.py b/python/snewpy/flavor_transformation.py index d30959ce..290273ea 100644 --- a/python/snewpy/flavor_transformation.py +++ b/python/snewpy/flavor_transformation.py @@ -1674,3 +1674,218 @@ def prob_xebar(self, t, E): else: return ( 1 - self.De3 - self.Ds3 ) / 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. From c089b09a069f3d112ccc08f07073b518ab28c7aa Mon Sep 17 00:00:00 2001 From: mvsantos001 Date: Fri, 14 Jul 2023 14:01:36 -0300 Subject: [PATCH 02/11] First version of QuantumDecoherence class --- python/snewpy/flavor_transformation.py | 121 +++++++++++++++++-------- 1 file changed, 82 insertions(+), 39 deletions(-) diff --git a/python/snewpy/flavor_transformation.py b/python/snewpy/flavor_transformation.py index 290273ea..bbc80f6f 100644 --- a/python/snewpy/flavor_transformation.py +++ b/python/snewpy/flavor_transformation.py @@ -1675,23 +1675,22 @@ def prob_xebar(self, t, E): return ( 1 - self.De3 - self.Ds3 ) / 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. +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, mass=1*u.eV/c.c**2, tau=1*u.day, dist=10*u.kpc, mh=MassHierarchy.NORMAL): + def __init__(self, mix_angles=None, Gamma3=1e-27*u.eV, Gamma8=1e-27*u.eV, 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. + 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. mh : MassHierarchy @@ -1712,26 +1711,79 @@ def __init__(self, mix_angles=None, mass=1*u.eV/c.c**2, tau=1*u.day, dist=10*u.k self.De2 = float((np.sin(theta12) * np.cos(theta13))**2) self.De3 = float(np.sin(theta13)**2) - self.m = mass - self.tau = tau + 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 + + def P11(self, E): + """Survival proability of state nu1 in vacuum. - def gamma(self, E): - """Decay width of the heaviest neutrino mass. + Parameters + ---------- + E : float + Energy. + + Returns + ------- + P11 : float + Survival proability of state nu1 in vacuum. + + :meta private: + """ + return 1/3 + 1/2 * np.exp(-(self.Gamma3 + self.Gamma8/3) * self.d) + 1/6 * np.exp(-self.Gamma8 * self.d) + + def P21(self, E): + """Transition proability from the state nu2 to nu1 in vacuum. Parameters ---------- E : float - Energy of the nu3. + Energy. Returns ------- - Gamma : float - Decay width of the neutrino mass, in units of 1/length. + P21 : float + Transition probability from the state nu2 to nu1 in vacuum. + Note that P21 = P12. :meta private: """ - return self.m*c.c / (E*self.tau) + return 1/3 - 1/2 * np.exp(-(self.Gamma3 + self.Gamma8/3) * self.d) + 1/6 * np.exp(-self.Gamma8 * self.d) + + def P31(self, E): + """Transition probability from the state nu3 to nu1 in vacuum. + + Parameters + ---------- + E : float + Energy. + + Returns + ------- + P31 : float + Transition proability from the state nu3 to nu1 in vacuum. + Note that P31 = P13. + + :meta private: + """ + return 1/3 - 1/3 * np.exp(-self.Gamma8 * self.d) + + def P33(self, E): + """Survival proability of state nu3 in vacuum. + + Parameters + ---------- + E : float + Energy. + + Returns + ------- + P33 : float + Survival proability of state nu3 in vacuum. + + :meta private: + """ + return 1/3 + 2/3 * np.exp(-self.Gamma8 * self.d) def prob_ee(self, t, E): """Electron neutrino survival probability. @@ -1750,12 +1802,10 @@ def prob_ee(self, t, E): """ # 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) + pe_array = self.P31(E)*self.De1 + self.P32(E)*self.De2 + self.P33(E)*self.De3 # IMO case. else: - pe_array = self.De2*np.exp(-self.gamma(E)*self.d) + \ - self.De3*(1-np.exp(-self.gamma(E)*self.d)) + 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): @@ -1773,12 +1823,7 @@ def prob_ex(self, t, E): 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 + return 1. - self.prob_ee(t,E) def prob_xx(self, t, E): """Flavor X neutrino survival probability. @@ -1829,7 +1874,13 @@ def prob_eebar(self, t, E): prob : float or ndarray Transition probability. """ - return self.De3 + # NMO case. + if self.mass_order == MassHierarchy.NORMAL: + pe_array = self.P11(E)*self.De1 + self.P12(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. @@ -1846,15 +1897,7 @@ def prob_exbar(self, t, E): 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 + return 1. - self.prob_eebar(t,E) def prob_xxbar(self, t, E): """X -> X antineutrino survival probability. @@ -1888,4 +1931,4 @@ def prob_xebar(self, t, E): prob : float or ndarray Transition probability. """ - return (1. - self.prob_eebar(t,E)) / 2. + return (1. - self.prob_eebar(t,E)) / 2. From 1ad741840e2bf1160307b862c139331dfdca22f5 Mon Sep 17 00:00:00 2001 From: mvsantos001 Date: Fri, 14 Jul 2023 15:03:36 -0300 Subject: [PATCH 03/11] Including P32 probability --- python/snewpy/flavor_transformation.py | 30 ++++++++++++++++++++------ 1 file changed, 24 insertions(+), 6 deletions(-) diff --git a/python/snewpy/flavor_transformation.py b/python/snewpy/flavor_transformation.py index bbc80f6f..03cc2eaf 100644 --- a/python/snewpy/flavor_transformation.py +++ b/python/snewpy/flavor_transformation.py @@ -1716,7 +1716,7 @@ def __init__(self, mix_angles=None, Gamma3=1e-27*u.eV, Gamma8=1e-27*u.eV, dist=1 self.d = dist def P11(self, E): - """Survival proability of state nu1 in vacuum. + """Survival probability of state nu1 in vacuum. Parameters ---------- @@ -1726,14 +1726,14 @@ def P11(self, E): Returns ------- P11 : float - Survival proability of state nu1 in vacuum. + Survival probability of state nu1 in vacuum. :meta private: """ return 1/3 + 1/2 * np.exp(-(self.Gamma3 + self.Gamma8/3) * self.d) + 1/6 * np.exp(-self.Gamma8 * self.d) def P21(self, E): - """Transition proability from the state nu2 to nu1 in vacuum. + """Transition probability from the state nu2 to nu1 in vacuum. Parameters ---------- @@ -1761,15 +1761,33 @@ def P31(self, E): Returns ------- P31 : float - Transition proability from the state nu3 to nu1 in vacuum. + 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 * 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 proability of state nu3 in vacuum. + """Survival probability of state nu3 in vacuum. Parameters ---------- @@ -1779,7 +1797,7 @@ def P33(self, E): Returns ------- P33 : float - Survival proability of state nu3 in vacuum. + Survival probability of state nu3 in vacuum. :meta private: """ From 8eddcb0ccfeb83b571841995a9d36c6013d016df Mon Sep 17 00:00:00 2001 From: mvsantos001 Date: Fri, 14 Jul 2023 15:11:32 -0300 Subject: [PATCH 04/11] Including P22 and P21 --- python/snewpy/flavor_transformation.py | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/python/snewpy/flavor_transformation.py b/python/snewpy/flavor_transformation.py index 03cc2eaf..9583354c 100644 --- a/python/snewpy/flavor_transformation.py +++ b/python/snewpy/flavor_transformation.py @@ -1750,6 +1750,24 @@ def P21(self, E): """ return 1/3 - 1/2 * np.exp(-(self.Gamma3 + self.Gamma8/3) * self.d) + 1/6 * np.exp(-self.Gamma8 * 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. @@ -1894,7 +1912,7 @@ def prob_eebar(self, t, E): """ # NMO case. if self.mass_order == MassHierarchy.NORMAL: - pe_array = self.P11(E)*self.De1 + self.P12(E)*self.De2 + self.P31(E)*self.De3 + 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 From 0151cc25b8cc364959f3b5996876cae674ea9a3c Mon Sep 17 00:00:00 2001 From: mvsantos001 Date: Fri, 14 Jul 2023 16:22:34 -0300 Subject: [PATCH 05/11] Including energy dependency on Gamma3 and Gamma8 --- python/snewpy/flavor_transformation.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/python/snewpy/flavor_transformation.py b/python/snewpy/flavor_transformation.py index 9583354c..c8d5c341 100644 --- a/python/snewpy/flavor_transformation.py +++ b/python/snewpy/flavor_transformation.py @@ -1680,7 +1680,7 @@ class QuantumDecoherence(FlavorTransformation): 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, mh=MassHierarchy.NORMAL): + 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 @@ -1693,6 +1693,11 @@ def __init__(self, mix_angles=None, Gamma3=1e-27*u.eV, Gamma8=1e-27*u.eV, dist=1 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. + E0 : astropy.units.quantity.Quantity + Reference energy in the power law Gamma = Gamma0*(E/E0)**n. mh : MassHierarchy MassHierarchy.NORMAL or MassHierarchy.INVERTED. """ @@ -1730,7 +1735,7 @@ def P11(self, E): :meta private: """ - return 1/3 + 1/2 * np.exp(-(self.Gamma3 + self.Gamma8/3) * self.d) + 1/6 * np.exp(-self.Gamma8 * self.d) + 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. @@ -1748,7 +1753,7 @@ def P21(self, E): :meta private: """ - return 1/3 - 1/2 * np.exp(-(self.Gamma3 + self.Gamma8/3) * self.d) + 1/6 * np.exp(-self.Gamma8 * self.d) + 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. @@ -1784,7 +1789,7 @@ def P31(self, E): :meta private: """ - return 1/3 - 1/3 * np.exp(-self.Gamma8 * self.d) + 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. @@ -1819,7 +1824,7 @@ def P33(self, E): :meta private: """ - return 1/3 + 2/3 * np.exp(-self.Gamma8 * self.d) + 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. From d1c93b98f9dfce91a2bf05edef816dfd07e3e92c Mon Sep 17 00:00:00 2001 From: mvsantos001 Date: Fri, 14 Jul 2023 16:27:02 -0300 Subject: [PATCH 06/11] Update flavor_transformation.py --- python/snewpy/flavor_transformation.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/python/snewpy/flavor_transformation.py b/python/snewpy/flavor_transformation.py index c8d5c341..908fdd7c 100644 --- a/python/snewpy/flavor_transformation.py +++ b/python/snewpy/flavor_transformation.py @@ -1719,6 +1719,8 @@ def __init__(self, mix_angles=None, Gamma3=1e-27*u.eV, Gamma8=1e-27*u.eV, dist=1 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. From 2211c015ecf8c302fa2564f67f74f1853a6b5665 Mon Sep 17 00:00:00 2001 From: mvsantos001 Date: Mon, 17 Jul 2023 08:36:43 -0300 Subject: [PATCH 07/11] Describe reference values of n and E0. --- python/snewpy/flavor_transformation.py | 6 ++++-- python/snewpy/test/test_04_xforms.py | 2 +- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/python/snewpy/flavor_transformation.py b/python/snewpy/flavor_transformation.py index 908fdd7c..b0f3667b 100644 --- a/python/snewpy/flavor_transformation.py +++ b/python/snewpy/flavor_transformation.py @@ -1695,9 +1695,11 @@ def __init__(self, mix_angles=None, Gamma3=1e-27*u.eV, Gamma8=1e-27*u.eV, dist=1 Distance to the supernova. n : float Exponent of power law for energy dependent quantum decoherence parameters, - i.e. Gamma = Gamma0*(E/E0)**n. + 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. + 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. """ diff --git a/python/snewpy/test/test_04_xforms.py b/python/snewpy/test/test_04_xforms.py index 676c0504..22c74673 100644 --- a/python/snewpy/test/test_04_xforms.py +++ b/python/snewpy/test/test_04_xforms.py @@ -7,7 +7,7 @@ AdiabaticMSW, NonAdiabaticMSWH, \ AdiabaticMSWes, NonAdiabaticMSWes, \ TwoFlavorDecoherence, ThreeFlavorDecoherence, \ - NeutrinoDecay + NeutrinoDecay, QuantumDecoherence from astropy import units as u from astropy import constants as c From 959a9dc2dbf87eaa15c509f72cf59364dd7e0753 Mon Sep 17 00:00:00 2001 From: mvsantos001 Date: Mon, 17 Jul 2023 08:40:15 -0300 Subject: [PATCH 08/11] Include QuantumDecoherence class in flavor_transformation_dict --- python/snewpy/snowglobes.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/snewpy/snowglobes.py b/python/snewpy/snowglobes.py index e70db700..50559098 100644 --- a/python/snewpy/snowglobes.py +++ b/python/snewpy/snowglobes.py @@ -73,7 +73,7 @@ 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_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), 'QuantumDecoherence_NMO': QuantumDecoherence(mh=MassHierarchy.NORMAL), 'QuantumDecoherence_IMO': QuantumDecoherence(mh=MassHierarchy.INVERTED)} flavor_transformation = flavor_transformation_dict[transformation_type] model_dir, model_file = os.path.split(os.path.abspath(model_path)) @@ -179,7 +179,7 @@ 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_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), 'QuantumDecoherence_NMO': QuantumDecoherence(mh=MassHierarchy.NORMAL), 'QuantumDecoherence_IMO': QuantumDecoherence(mh=MassHierarchy.INVERTED)} flavor_transformation = flavor_transformation_dict[transformation_type] model_dir, model_file = os.path.split(os.path.abspath(model_path)) From 67b79eb611802887d6039dd81784991f3dafab69 Mon Sep 17 00:00:00 2001 From: mvsantos001 Date: Mon, 17 Jul 2023 19:10:09 -0300 Subject: [PATCH 09/11] Implementation of qd in test_04... --- python/snewpy/test/test_04_xforms.py | 130 +++++++++++++++++++++++++++ 1 file changed, 130 insertions(+) diff --git a/python/snewpy/test/test_04_xforms.py b/python/snewpy/test/test_04_xforms.py index 22c74673..360cee7d 100644 --- a/python/snewpy/test/test_04_xforms.py +++ b/python/snewpy/test/test_04_xforms.py @@ -32,6 +32,12 @@ def setUp(self): self.lifetime = 1 * u.day self.distance = 10 * u.kpc + # Quantum decoherence parameters; see arXiv:2306.17591. + self.gamma3 = 1e-27 + self.gamma8 = 1e-27 + self.n = 0 + self.energy_ref = 10 * u.MeV + def test_noxform(self): """ Survival probabilities for no oscillations @@ -486,3 +492,127 @@ def test_nudecay_imo_default_mixing(self): self.assertTrue(np.array_equal(xform.prob_exbar(self.t, self.E), prob_exbar)) self.assertTrue(np.array_equal(xform.prob_xxbar(self.t, self.E), 1. - 0.5*prob_exbar)) self.assertEqual(xform.prob_xebar(self.t, self.E), 0.5*(1. - De3)) + + def test_quantumdecoherence_nmo(self): + """ + Quantum Decoherence with NMO and new mixing angles + """ + # Override the default mixing angles. + xform = QuantumDecoherence(mix_angles=(self.theta12, self.theta13, self.theta23), Gamma3=self.gamma3, Gamma8=self.gamma8, dist=self.distance, n=self.n, E0=self.energy_ref, mh=MassHierarchy.NORMAL) + + # Test computation survival and transition probabilities of mass states. + _E = 10*u.MeV + self.assertTrue(xform.P11(_E) == 1/3 + 1/2 * np.exp(-(self.gamma3 * (_E/self.energy_ref)**self.n + self.gamma8 * (_E/self.energy_ref)**self.n / 3) * self.distance) + 1/6 * np.exp(-self.gamma8 * (_E/self.energy_ref)**self.n * self.distance)) + self.assertTrue(xform.P21(_E) == 1/3 - 1/2 * np.exp(-(self.gamma3 * (_E/self.energy_ref)**self.n + self.gamma8 * (_E/self.energy_ref)**self.n / 3) * self.distance) + 1/6 * np.exp(-self.gamma8 * (_E/self.energy_ref)**self.n * self.distance)) + self.assertTrue(xform.P22(_E) == 1/3 + 1/2 * np.exp(-(self.gamma3 * (_E/self.energy_ref)**self.n + self.gamma8 * (_E/self.energy_ref)**self.n / 3) * self.distance) + 1/6 * np.exp(-self.gamma8 * (_E/self.energy_ref)**self.n * self.distance)) + self.assertTrue(xform.P31(_E) == 1/3 - 1/3 * np.exp(-self.gamma8 * (_E/self.energy_ref)**self.n * self.distance)) + self.assertTrue(xform.P32(_E) == 1/3 - 1/3 * np.exp(-self.gamma8 * (_E/self.energy_ref)**self.n * self.distance)) + self.assertTrue(xform.P33(_E) == 1/3 + 2/3 * np.exp(-self.gamma8 * (_E/self.E0)**self.n * self.distance)) + + De1 = (cos(self.theta12) * cos(self.theta13))**2 + De2 = (sin(self.theta12) * cos(self.theta13))**2 + De3 = sin(self.theta13)**2 + + # Check flavor transition probabilities. + prob_ee = np.asarray([xform.P31(_E)*De1 + xform.P32(_E)*De2 + xform.P33(_E)*De3] for _E in self.E) + + self.assertTrue(np.array_equal(xform.prob_ee(self.t, self.E), prob_ee)) + self.assertEqual(xform.prob_ex(self.t, self.E), 1 - prob_ee) + self.assertEqual(xform.prob_xx(self.t, self.E), 1 - 0.5*(1 - prob_ee)) + self.assertTrue(np.array_equal(xform.prob_xe(self.t, self.E), 0.5*(1 - prob_ee))) + + prob_eebar = np.asarray([xform.P11(_E)*De1 + xform.P21(_E)*De2 + xform.P31(_E)*De3 for _E in self.E]) + + self.assertEqual(xform.prob_exbar(self.t, self.E), 1 - prob_eebar) + self.assertTrue(np.array_equal(xform.prob_xxbar(self.t, self.E), 1. - 0.5*(1 - prob_eebar))) + self.assertEqual(xform.prob_xebar(self.t, self.E), 0.5*(1. - prob_eebar)) + + def test_quantumdecoherence_nmo_default_mixing(self): + """ + Quantum decoherence with NMO and default mixing angles + """ + # Use default mixing angles defined in the submodule. + xform = QuantumDecoherence(Gamma3=self.gamma3, Gamma8=self.gamma8, dist=self.distance, n=self.n, E0=self.energy_ref) + + # Check transition probabilities (normal hierarchy is default). + mixpars = MixingParameters() + th12, th13, th23 = mixpars.get_mixing_angles() + + De1 = (cos(th12) * cos(th13))**2 + De2 = (sin(th12) * cos(th13))**2 + De3 = sin(th13)**2 + + prob_ee = np.asarray([xform.P31(_E)*De1 + xform.P32(_E)*De2 + xform.P33(_E)*De3] for _E in self.E) + + self.assertTrue(np.array_equal(xform.prob_ee(self.t, self.E), prob_ee)) + self.assertEqual(xform.prob_ex(self.t, self.E), 1 - prob_ee) + self.assertEqual(xform.prob_xx(self.t, self.E), 1 - 0.5*(1 - prob_ee)) + self.assertTrue(np.array_equal(xform.prob_xe(self.t, self.E), 0.5*(1 - prob_ee))) + + prob_eebar = np.asarray([xform.P11(_E)*De1 + xform.P21(_E)*De2 + xform.P31(_E)*De3 for _E in self.E]) + + self.assertEqual(xform.prob_exbar(self.t, self.E), 1 - prob_eebar) + self.assertTrue(np.array_equal(xform.prob_xxbar(self.t, self.E), 1. - 0.5*(1 - prob_eebar))) + self.assertEqual(xform.prob_xebar(self.t, self.E), 0.5*(1. - prob_eebar)) + + def test_quantumdecoherence_imo(self): + """ + Quantum Decoherence with IMO and new mixing angles + """ + # Override the default mixing angles. + xform = QuantumDecoherence(mix_angles=(self.theta12, self.theta13, self.theta23), Gamma3=self.gamma3, Gamma8=self.gamma8, dist=self.distance, n=self.n, E0=self.energy_ref, mh=MassHierarchy.INVERTED) + + # Test computation survival and transition probabilities of mass states. + _E = 10*u.MeV + self.assertTrue(xform.P11(_E) == 1/3 + 1/2 * np.exp(-(self.gamma3 * (_E/self.energy_ref)**self.n + self.gamma8 * (_E/self.energy_ref)**self.n / 3) * self.distance) + 1/6 * np.exp(-self.gamma8 * (_E/self.energy_ref)**self.n * self.distance)) + self.assertTrue(xform.P21(_E) == 1/3 - 1/2 * np.exp(-(self.gamma3 * (_E/self.energy_ref)**self.n + self.gamma8 * (_E/self.energy_ref)**self.n / 3) * self.distance) + 1/6 * np.exp(-self.gamma8 * (_E/self.energy_ref)**self.n * self.distance)) + self.assertTrue(xform.P22(_E) == 1/3 + 1/2 * np.exp(-(self.gamma3 * (_E/self.energy_ref)**self.n + self.gamma8 * (_E/self.energy_ref)**self.n / 3) * self.distance) + 1/6 * np.exp(-self.gamma8 * (_E/self.energy_ref)**self.n * self.distance)) + self.assertTrue(xform.P31(_E) == 1/3 - 1/3 * np.exp(-self.gamma8 * (_E/self.energy_ref)**self.n * self.distance)) + self.assertTrue(xform.P32(_E) == 1/3 - 1/3 * np.exp(-self.gamma8 * (_E/self.energy_ref)**self.n * self.distance)) + self.assertTrue(xform.P33(_E) == 1/3 + 2/3 * np.exp(-self.gamma8 * (_E/self.E0)**self.n * self.distance)) + + De1 = (cos(self.theta12) * cos(self.theta13))**2 + De2 = (sin(self.theta12) * cos(self.theta13))**2 + De3 = sin(self.theta13)**2 + + # Check transition probabilities. + prob_ee = np.asarray([xform.P22(_E)*De2 + xform.P21(_E)*De1 + xform.P32(_E)*De3] for _E in self.E) + + self.assertTrue(np.array_equal(xform.prob_ee(self.t, self.E), prob_ee)) + self.assertEqual(xform.prob_ex(self.t, self.E), 1 - prob_ee) + self.assertEqual(xform.prob_xx(self.t, self.E), 1 - 0.5*(1 - prob_ee)) + self.assertTrue(np.array_equal(xform.prob_xe(self.t, self.E), 0.5*(1 - prob_ee))) + + prob_eebar = np.asarray([xform.P31(_E)*De1 + xform.P32(_E)*De2 + xform.P33(_E)*De3 for _E in self.E]) + + self.assertEqual(xform.prob_exbar(self.t, self.E), 1 - prob_eebar) + self.assertTrue(np.array_equal(xform.prob_xxbar(self.t, self.E), 1. - 0.5*(1 - prob_eebar))) + self.assertEqual(xform.prob_xebar(self.t, self.E), 0.5*(1. - prob_eebar)) + + def test_quantumdecoherence_imo_default_mixing(self): + """ + Quantum decoherence with IMO and default mixing angles + """ + # Use default mixing angles defined in the submodule. + xform = QuantumDecoherence(Gamma3=self.gamma3, Gamma8=self.gamma8, dist=self.distance, n=self.n, E0=self.energy_ref, mh=MassHierarchy.INVERTED) + + # Check transition probabilities. + mixpars = MixingParameters() + th12, th13, th23 = mixpars.get_mixing_angles() + + De1 = (cos(th12) * cos(th13))**2 + De2 = (sin(th12) * cos(th13))**2 + De3 = sin(th13)**2 + + prob_ee = np.asarray([xform.P22(_E)*De2 + xform.P21(_E)*De1 + xform.P32(_E)*De3] for _E in self.E) + + self.assertTrue(np.array_equal(xform.prob_ee(self.t, self.E), prob_ee)) + self.assertEqual(xform.prob_ex(self.t, self.E), 1 - prob_ee) + self.assertEqual(xform.prob_xx(self.t, self.E), 1 - 0.5*(1 - prob_ee)) + self.assertTrue(np.array_equal(xform.prob_xe(self.t, self.E), 0.5*(1 - prob_ee))) + + prob_eebar = np.asarray([xform.P31(_E)*De1 + xform.P32(_E)*De2 + xform.P33(_E)*De3 for _E in self.E]) + + self.assertEqual(xform.prob_exbar(self.t, self.E), 1 - prob_eebar) + self.assertTrue(np.array_equal(xform.prob_xxbar(self.t, self.E), 1. - 0.5*(1 - prob_eebar))) + self.assertEqual(xform.prob_xebar(self.t, self.E), 0.5*(1. - prob_eebar)) \ No newline at end of file From 4f6e8f9132e4e69392513b7c67a8e1a9d2d5325e Mon Sep 17 00:00:00 2001 From: mvsantos001 Date: Thu, 19 Oct 2023 09:49:38 -0300 Subject: [PATCH 10/11] Solving an issue with units in test_04_xforms.py --- python/snewpy/test/test_04_xforms.py | 58 ++++++++++++++-------------- 1 file changed, 29 insertions(+), 29 deletions(-) diff --git a/python/snewpy/test/test_04_xforms.py b/python/snewpy/test/test_04_xforms.py index 360cee7d..365fbbb8 100644 --- a/python/snewpy/test/test_04_xforms.py +++ b/python/snewpy/test/test_04_xforms.py @@ -33,8 +33,8 @@ def setUp(self): self.distance = 10 * u.kpc # Quantum decoherence parameters; see arXiv:2306.17591. - self.gamma3 = 1e-27 - self.gamma8 = 1e-27 + self.gamma3 = (1e-27 * u.eV / (c.hbar.to('eV s') * c.c)).to('1/kpc') + self.gamma8 = (1e-27 * u.eV / (c.hbar.to('eV s') * c.c)).to('1/kpc') self.n = 0 self.energy_ref = 10 * u.MeV @@ -498,7 +498,7 @@ def test_quantumdecoherence_nmo(self): Quantum Decoherence with NMO and new mixing angles """ # Override the default mixing angles. - xform = QuantumDecoherence(mix_angles=(self.theta12, self.theta13, self.theta23), Gamma3=self.gamma3, Gamma8=self.gamma8, dist=self.distance, n=self.n, E0=self.energy_ref, mh=MassHierarchy.NORMAL) + xform = QuantumDecoherence(mix_angles=(self.theta12, self.theta13, self.theta23), Gamma3=self.gamma3 * (c.hbar.to('eV s') * c.c.to('kpc / s')), Gamma8=self.gamma8 * (c.hbar.to('eV s') * c.c.to('kpc / s')), dist=self.distance, n=self.n, E0=self.energy_ref, mh=MassHierarchy.NORMAL) # Test computation survival and transition probabilities of mass states. _E = 10*u.MeV @@ -507,32 +507,32 @@ def test_quantumdecoherence_nmo(self): self.assertTrue(xform.P22(_E) == 1/3 + 1/2 * np.exp(-(self.gamma3 * (_E/self.energy_ref)**self.n + self.gamma8 * (_E/self.energy_ref)**self.n / 3) * self.distance) + 1/6 * np.exp(-self.gamma8 * (_E/self.energy_ref)**self.n * self.distance)) self.assertTrue(xform.P31(_E) == 1/3 - 1/3 * np.exp(-self.gamma8 * (_E/self.energy_ref)**self.n * self.distance)) self.assertTrue(xform.P32(_E) == 1/3 - 1/3 * np.exp(-self.gamma8 * (_E/self.energy_ref)**self.n * self.distance)) - self.assertTrue(xform.P33(_E) == 1/3 + 2/3 * np.exp(-self.gamma8 * (_E/self.E0)**self.n * self.distance)) + self.assertAlmostEqual(float(xform.P33(_E)), float(1/3 + 2/3 * np.exp(-self.gamma8 * (_E/self.energy_ref)**self.n * self.distance)), places=12) De1 = (cos(self.theta12) * cos(self.theta13))**2 De2 = (sin(self.theta12) * cos(self.theta13))**2 De3 = sin(self.theta13)**2 # Check flavor transition probabilities. - prob_ee = np.asarray([xform.P31(_E)*De1 + xform.P32(_E)*De2 + xform.P33(_E)*De3] for _E in self.E) + prob_ee = np.asarray([xform.P31(_E)*De1 + xform.P32(_E)*De2 + xform.P33(_E)*De3 for _E in self.E]) self.assertTrue(np.array_equal(xform.prob_ee(self.t, self.E), prob_ee)) - self.assertEqual(xform.prob_ex(self.t, self.E), 1 - prob_ee) - self.assertEqual(xform.prob_xx(self.t, self.E), 1 - 0.5*(1 - prob_ee)) + self.assertTrue(np.array_equal(xform.prob_ex(self.t, self.E), 1 - prob_ee)) + self.assertTrue(np.array_equal(xform.prob_xx(self.t, self.E), 1 - 0.5*(1 - prob_ee))) self.assertTrue(np.array_equal(xform.prob_xe(self.t, self.E), 0.5*(1 - prob_ee))) prob_eebar = np.asarray([xform.P11(_E)*De1 + xform.P21(_E)*De2 + xform.P31(_E)*De3 for _E in self.E]) - self.assertEqual(xform.prob_exbar(self.t, self.E), 1 - prob_eebar) + self.assertTrue(np.array_equal(xform.prob_exbar(self.t, self.E), 1 - prob_eebar)) self.assertTrue(np.array_equal(xform.prob_xxbar(self.t, self.E), 1. - 0.5*(1 - prob_eebar))) - self.assertEqual(xform.prob_xebar(self.t, self.E), 0.5*(1. - prob_eebar)) + self.assertTrue(np.array_equal(xform.prob_xebar(self.t, self.E), 0.5*(1. - prob_eebar))) def test_quantumdecoherence_nmo_default_mixing(self): """ Quantum decoherence with NMO and default mixing angles """ # Use default mixing angles defined in the submodule. - xform = QuantumDecoherence(Gamma3=self.gamma3, Gamma8=self.gamma8, dist=self.distance, n=self.n, E0=self.energy_ref) + xform = QuantumDecoherence(Gamma3=self.gamma3 * (c.hbar.to('eV s') * c.c.to('kpc / s')), Gamma8=self.gamma8 * (c.hbar.to('eV s') * c.c.to('kpc / s')), dist=self.distance, n=self.n, E0=self.energy_ref) # Check transition probabilities (normal hierarchy is default). mixpars = MixingParameters() @@ -542,25 +542,25 @@ def test_quantumdecoherence_nmo_default_mixing(self): De2 = (sin(th12) * cos(th13))**2 De3 = sin(th13)**2 - prob_ee = np.asarray([xform.P31(_E)*De1 + xform.P32(_E)*De2 + xform.P33(_E)*De3] for _E in self.E) + prob_ee = np.asarray([xform.P31(_E)*De1 + xform.P32(_E)*De2 + xform.P33(_E)*De3 for _E in self.E]) self.assertTrue(np.array_equal(xform.prob_ee(self.t, self.E), prob_ee)) - self.assertEqual(xform.prob_ex(self.t, self.E), 1 - prob_ee) - self.assertEqual(xform.prob_xx(self.t, self.E), 1 - 0.5*(1 - prob_ee)) + self.assertTrue(np.array_equal(xform.prob_ex(self.t, self.E), 1 - prob_ee)) + self.assertTrue(np.array_equal(xform.prob_xx(self.t, self.E), 1 - 0.5*(1 - prob_ee))) self.assertTrue(np.array_equal(xform.prob_xe(self.t, self.E), 0.5*(1 - prob_ee))) prob_eebar = np.asarray([xform.P11(_E)*De1 + xform.P21(_E)*De2 + xform.P31(_E)*De3 for _E in self.E]) - self.assertEqual(xform.prob_exbar(self.t, self.E), 1 - prob_eebar) + self.assertTrue(np.array_equal(xform.prob_exbar(self.t, self.E), 1 - prob_eebar)) self.assertTrue(np.array_equal(xform.prob_xxbar(self.t, self.E), 1. - 0.5*(1 - prob_eebar))) - self.assertEqual(xform.prob_xebar(self.t, self.E), 0.5*(1. - prob_eebar)) + self.assertTrue(np.array_equal(xform.prob_xebar(self.t, self.E), 0.5*(1. - prob_eebar))) def test_quantumdecoherence_imo(self): """ Quantum Decoherence with IMO and new mixing angles """ # Override the default mixing angles. - xform = QuantumDecoherence(mix_angles=(self.theta12, self.theta13, self.theta23), Gamma3=self.gamma3, Gamma8=self.gamma8, dist=self.distance, n=self.n, E0=self.energy_ref, mh=MassHierarchy.INVERTED) + xform = QuantumDecoherence(mix_angles=(self.theta12, self.theta13, self.theta23), Gamma3=self.gamma3 * (c.hbar.to('eV s') * c.c.to('kpc / s')), Gamma8=self.gamma8 * (c.hbar.to('eV s') * c.c.to('kpc / s')), dist=self.distance, n=self.n, E0=self.energy_ref, mh=MassHierarchy.INVERTED) # Test computation survival and transition probabilities of mass states. _E = 10*u.MeV @@ -569,50 +569,50 @@ def test_quantumdecoherence_imo(self): self.assertTrue(xform.P22(_E) == 1/3 + 1/2 * np.exp(-(self.gamma3 * (_E/self.energy_ref)**self.n + self.gamma8 * (_E/self.energy_ref)**self.n / 3) * self.distance) + 1/6 * np.exp(-self.gamma8 * (_E/self.energy_ref)**self.n * self.distance)) self.assertTrue(xform.P31(_E) == 1/3 - 1/3 * np.exp(-self.gamma8 * (_E/self.energy_ref)**self.n * self.distance)) self.assertTrue(xform.P32(_E) == 1/3 - 1/3 * np.exp(-self.gamma8 * (_E/self.energy_ref)**self.n * self.distance)) - self.assertTrue(xform.P33(_E) == 1/3 + 2/3 * np.exp(-self.gamma8 * (_E/self.E0)**self.n * self.distance)) + self.assertAlmostEqual(float(xform.P33(_E)), float(1/3 + 2/3 * np.exp(-self.gamma8 * (_E/self.energy_ref)**self.n * self.distance)), places=12) De1 = (cos(self.theta12) * cos(self.theta13))**2 De2 = (sin(self.theta12) * cos(self.theta13))**2 De3 = sin(self.theta13)**2 # Check transition probabilities. - prob_ee = np.asarray([xform.P22(_E)*De2 + xform.P21(_E)*De1 + xform.P32(_E)*De3] for _E in self.E) + prob_ee = np.asarray([xform.P22(_E)*De2 + xform.P21(_E)*De1 + xform.P32(_E)*De3 for _E in self.E]) self.assertTrue(np.array_equal(xform.prob_ee(self.t, self.E), prob_ee)) - self.assertEqual(xform.prob_ex(self.t, self.E), 1 - prob_ee) - self.assertEqual(xform.prob_xx(self.t, self.E), 1 - 0.5*(1 - prob_ee)) + self.assertTrue(np.array_equal(xform.prob_ex(self.t, self.E), 1 - prob_ee)) + self.assertTrue(np.array_equal(xform.prob_xx(self.t, self.E), 1 - 0.5*(1 - prob_ee))) self.assertTrue(np.array_equal(xform.prob_xe(self.t, self.E), 0.5*(1 - prob_ee))) prob_eebar = np.asarray([xform.P31(_E)*De1 + xform.P32(_E)*De2 + xform.P33(_E)*De3 for _E in self.E]) - self.assertEqual(xform.prob_exbar(self.t, self.E), 1 - prob_eebar) + self.assertTrue(np.array_equal(xform.prob_exbar(self.t, self.E), 1 - prob_eebar)) self.assertTrue(np.array_equal(xform.prob_xxbar(self.t, self.E), 1. - 0.5*(1 - prob_eebar))) - self.assertEqual(xform.prob_xebar(self.t, self.E), 0.5*(1. - prob_eebar)) + self.assertTrue(np.array_equal(xform.prob_xebar(self.t, self.E), 0.5*(1. - prob_eebar))) def test_quantumdecoherence_imo_default_mixing(self): """ Quantum decoherence with IMO and default mixing angles """ # Use default mixing angles defined in the submodule. - xform = QuantumDecoherence(Gamma3=self.gamma3, Gamma8=self.gamma8, dist=self.distance, n=self.n, E0=self.energy_ref, mh=MassHierarchy.INVERTED) + xform = QuantumDecoherence(Gamma3=self.gamma3 * (c.hbar.to('eV s') * c.c.to('kpc / s')), Gamma8=self.gamma8 * (c.hbar.to('eV s') * c.c.to('kpc / s')), dist=self.distance, n=self.n, E0=self.energy_ref, mh=MassHierarchy.INVERTED) # Check transition probabilities. - mixpars = MixingParameters() + mixpars = MixingParameters(MassHierarchy.INVERTED) th12, th13, th23 = mixpars.get_mixing_angles() De1 = (cos(th12) * cos(th13))**2 De2 = (sin(th12) * cos(th13))**2 De3 = sin(th13)**2 - prob_ee = np.asarray([xform.P22(_E)*De2 + xform.P21(_E)*De1 + xform.P32(_E)*De3] for _E in self.E) + prob_ee = np.asarray([xform.P22(_E)*De2 + xform.P21(_E)*De1 + xform.P32(_E)*De3 for _E in self.E]) self.assertTrue(np.array_equal(xform.prob_ee(self.t, self.E), prob_ee)) - self.assertEqual(xform.prob_ex(self.t, self.E), 1 - prob_ee) - self.assertEqual(xform.prob_xx(self.t, self.E), 1 - 0.5*(1 - prob_ee)) + self.assertTrue(np.array_equal(xform.prob_ex(self.t, self.E), 1 - prob_ee)) + self.assertTrue(np.array_equal(xform.prob_xx(self.t, self.E), 1 - 0.5*(1 - prob_ee))) self.assertTrue(np.array_equal(xform.prob_xe(self.t, self.E), 0.5*(1 - prob_ee))) prob_eebar = np.asarray([xform.P31(_E)*De1 + xform.P32(_E)*De2 + xform.P33(_E)*De3 for _E in self.E]) - self.assertEqual(xform.prob_exbar(self.t, self.E), 1 - prob_eebar) + self.assertTrue(np.array_equal(xform.prob_exbar(self.t, self.E), 1 - prob_eebar)) self.assertTrue(np.array_equal(xform.prob_xxbar(self.t, self.E), 1. - 0.5*(1 - prob_eebar))) - self.assertEqual(xform.prob_xebar(self.t, self.E), 0.5*(1. - prob_eebar)) \ No newline at end of file + self.assertTrue(np.array_equal(xform.prob_xebar(self.t, self.E), 0.5*(1. - prob_eebar))) \ No newline at end of file From 52fa478dadca5824d138bdd5f49291462e8bd7c7 Mon Sep 17 00:00:00 2001 From: Jost Migenda Date: Wed, 22 Nov 2023 23:34:02 +0000 Subject: [PATCH 11/11] Clean up some unnecessary unit conversions --- python/snewpy/test/test_04_xforms.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/python/snewpy/test/test_04_xforms.py b/python/snewpy/test/test_04_xforms.py index 365fbbb8..0ed0f99e 100644 --- a/python/snewpy/test/test_04_xforms.py +++ b/python/snewpy/test/test_04_xforms.py @@ -33,8 +33,8 @@ def setUp(self): self.distance = 10 * u.kpc # Quantum decoherence parameters; see arXiv:2306.17591. - self.gamma3 = (1e-27 * u.eV / (c.hbar.to('eV s') * c.c)).to('1/kpc') - self.gamma8 = (1e-27 * u.eV / (c.hbar.to('eV s') * c.c)).to('1/kpc') + 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 @@ -498,7 +498,7 @@ def test_quantumdecoherence_nmo(self): Quantum Decoherence with NMO and new mixing angles """ # Override the default mixing angles. - xform = QuantumDecoherence(mix_angles=(self.theta12, self.theta13, self.theta23), Gamma3=self.gamma3 * (c.hbar.to('eV s') * c.c.to('kpc / s')), Gamma8=self.gamma8 * (c.hbar.to('eV s') * c.c.to('kpc / s')), dist=self.distance, n=self.n, E0=self.energy_ref, mh=MassHierarchy.NORMAL) + xform = QuantumDecoherence(mix_angles=(self.theta12, self.theta13, self.theta23), Gamma3=self.gamma3 * c.hbar * c.c, Gamma8=self.gamma8 * c.hbar * c.c, dist=self.distance, n=self.n, E0=self.energy_ref, mh=MassHierarchy.NORMAL) # Test computation survival and transition probabilities of mass states. _E = 10*u.MeV @@ -532,7 +532,7 @@ def test_quantumdecoherence_nmo_default_mixing(self): Quantum decoherence with NMO and default mixing angles """ # Use default mixing angles defined in the submodule. - xform = QuantumDecoherence(Gamma3=self.gamma3 * (c.hbar.to('eV s') * c.c.to('kpc / s')), Gamma8=self.gamma8 * (c.hbar.to('eV s') * c.c.to('kpc / s')), dist=self.distance, n=self.n, E0=self.energy_ref) + xform = QuantumDecoherence(Gamma3=self.gamma3 * c.hbar * c.c, Gamma8=self.gamma8 * c.hbar * c.c, dist=self.distance, n=self.n, E0=self.energy_ref) # Check transition probabilities (normal hierarchy is default). mixpars = MixingParameters() @@ -560,7 +560,7 @@ def test_quantumdecoherence_imo(self): Quantum Decoherence with IMO and new mixing angles """ # Override the default mixing angles. - xform = QuantumDecoherence(mix_angles=(self.theta12, self.theta13, self.theta23), Gamma3=self.gamma3 * (c.hbar.to('eV s') * c.c.to('kpc / s')), Gamma8=self.gamma8 * (c.hbar.to('eV s') * c.c.to('kpc / s')), dist=self.distance, n=self.n, E0=self.energy_ref, mh=MassHierarchy.INVERTED) + xform = QuantumDecoherence(mix_angles=(self.theta12, self.theta13, self.theta23), Gamma3=self.gamma3 * c.hbar * c.c, Gamma8=self.gamma8 * c.hbar * c.c, dist=self.distance, n=self.n, E0=self.energy_ref, mh=MassHierarchy.INVERTED) # Test computation survival and transition probabilities of mass states. _E = 10*u.MeV @@ -594,7 +594,7 @@ def test_quantumdecoherence_imo_default_mixing(self): Quantum decoherence with IMO and default mixing angles """ # Use default mixing angles defined in the submodule. - xform = QuantumDecoherence(Gamma3=self.gamma3 * (c.hbar.to('eV s') * c.c.to('kpc / s')), Gamma8=self.gamma8 * (c.hbar.to('eV s') * c.c.to('kpc / s')), dist=self.distance, n=self.n, E0=self.energy_ref, mh=MassHierarchy.INVERTED) + xform = QuantumDecoherence(Gamma3=self.gamma3 * c.hbar * c.c, Gamma8=self.gamma8 * c.hbar * c.c, dist=self.distance, n=self.n, E0=self.energy_ref, mh=MassHierarchy.INVERTED) # Check transition probabilities. mixpars = MixingParameters(MassHierarchy.INVERTED)