From 4675877adf45edbe2b5065ef903e540a5407de1f Mon Sep 17 00:00:00 2001 From: T Ehrhardt Date: Fri, 7 Nov 2025 12:30:27 +0100 Subject: [PATCH 1/7] introduce option to use NLO CC matter potential in prob3 + superficial mods: lazy logging, global overview of physics parameterisation choices, catch case of tomography_type string not being among accepted choices --- pisa/stages/osc/prob3.py | 120 ++++++++++++++++++++++++++++----------- 1 file changed, 86 insertions(+), 34 deletions(-) diff --git a/pisa/stages/osc/prob3.py b/pisa/stages/osc/prob3.py index 4f9023a14..7518f0d38 100644 --- a/pisa/stages/osc/prob3.py +++ b/pisa/stages/osc/prob3.py @@ -23,15 +23,52 @@ from pisa.stages.osc.prob3numba.numba_osc_hostfuncs import propagate_array, fill_probs from pisa.utils.resources import find_resource -__all__ = ['prob3', 'init_test'] +__all__ = ['prob3', 'init_test', + 'LRI_TYPES', 'NSI_TYPES', 'TOMOGRAPHY_TYPES'] + +LRI_TYPES = ['emu-symmetry', 'etau-symmetry', 'mutau-symmetry'] + +NSI_TYPES = ['standard', 'vacuum-like'] + +TOMOGRAPHY_TYPES = ['mass_of_earth', 'mass_of_core_w_constrain', 'mass_of_core_wo_constrain'] class prob3(Stage): # pylint: disable=invalid-name """ - Prob3-like oscillation PISA Pi class + Extended Prob3-like oscillations class Parameters ---------- + + include_nlo : bool (default: `False`) + Whether to include a +2.0% NLO correction to the SM CC + matter potential, as per + https://inspirehep.net/literature/2914951. + + nsi_type : str or `None` (default: `None`) + Choice of propagation/NC NSI parameterization. + If string, either 'standard' or 'vacuum-like'. + + reparam_mix_matrix : bool (default: `False`) + Whether to rephase the parameterization of the + leptonic mixing matrix from its PDG default by + diag(e^(i*delta_CP), 1, 1). This has no impact on + oscillation probabilities in the *absence* of NSI. + TODO: motivation + + neutrino_decay : bool (default: `False`) + Whether to invoke neutrino decay with oscillations. + TODO: details + + tomography_type : str or `None` (default: `None`) + If string, either 'mass_of_earth', 'mass_of_core_w_constrain', + or 'mass_of_core_wo_constrain'. + + lri_type : str or `None` (default: `None`) + Choice of model/parameterisation of long-range interactions. + If string, either 'emu-symmetry', 'etau-symmetry', + or 'mutau-symmetry'. + params Expected params .. :: @@ -88,6 +125,7 @@ class prob3(Stage): # pylint: disable=invalid-name def __init__( self, + include_nlo=False, nsi_type=None, reparam_mix_matrix=False, neutrino_decay=False, @@ -120,10 +158,13 @@ def __init__( 'weights' ) + self.include_nlo = include_nlo + """Whether to include a +2.0% NLO correction to the + SM CC matter potential.""" # Check whether and if so with which NSI parameters we are to work. if nsi_type is not None: - choices = ['standard', 'vacuum-like'] + choices = NSI_TYPES nsi_type = nsi_type.strip().lower() if not nsi_type in choices: raise ValueError( @@ -132,7 +173,7 @@ def __init__( ) self.nsi_type = nsi_type """Type of NSI to assume.""" - self.tomography_type = tomography_type + self.reparam_mix_matrix = reparam_mix_matrix """Use a PMNS mixing matrix parameterisation that differs from the standard one by an overall phase matrix @@ -178,7 +219,7 @@ def __init__( decay_params = () if lri_type is not None: - choices = ['emu-symmetry', 'etau-symmetry', 'mutau-symmetry'] + choices = LRI_TYPES lri_type = lri_type.strip().lower() if not lri_type in choices: raise ValueError( @@ -193,17 +234,26 @@ def __init__( lri_params = ('v_lri',) - if self.tomography_type is None: + if tomography_type is None: tomography_params = () - elif self.tomography_type == 'mass_of_earth': - tomography_params = ('density_scale',) - elif self.tomography_type == 'mass_of_core_w_constrain': - tomography_params = ('core_density_scale',) - elif self.tomography_type == 'mass_of_core_wo_constrain': - tomography_params = ('core_density_scale', - 'innermantle_density_scale', - 'middlemantle_density_scale' - ) + else: + tomography_type = tomography_type.strip().lower() + choices = TOMOGRAPHY_TYPES + if not tomography_type in choices: + raise ValueError( + 'Chosen tomography type "%s" not available! Choose one of %s.' + % (tomography_type, choices) + ) + if tomography_type == 'mass_of_earth': + tomography_params = ('density_scale',) + elif tomography_type == 'mass_of_core_w_constrain': + tomography_params = ('core_density_scale',) + elif tomography_type == 'mass_of_core_wo_constrain': + tomography_params = ('core_density_scale', + 'innermantle_density_scale', + 'middlemantle_density_scale' + ) + self.tomography_type = tomography_type expected_params = (expected_params + nsi_params + decay_params @@ -216,7 +266,6 @@ def __init__( **std_kwargs, ) - self.layers = None self.osc_params = None self.nsi_params = None @@ -231,9 +280,9 @@ def __init__( # (NSI), which is why we can simply treat it as a constant here. self.gen_mat_pot_matrix_complex = None """Interaction Hamiltonian without the factor sqrt(2)*G_F*N_e.""" - self.YeI = None - self.YeO = None - self.YeM = None + self.YeI = None # pylint: disable=invalid-name + self.YeO = None # pylint: disable=invalid-name + self.YeM = None # pylint: disable=invalid-name def setup_function(self): @@ -241,7 +290,7 @@ def setup_function(self): self.osc_params = OscParams() if self.reparam_mix_matrix: logging.debug( - 'Working with reparameterizated version of mixing matrix.' + 'Working with reparameterized version of mixing matrix.' ) else: logging.debug( @@ -265,13 +314,13 @@ def setup_function(self): if self.tomography_type == "mass_of_earth": - logging.debug('Working with a single density scaling factor.') + logging.debug('Working with tomography with a single density scaling factor.') self.tomography_params = Mass_scaling() elif self.tomography_type == "mass_of_core_w_constrain": - logging.debug('Working with different scaling for different layers.') + logging.debug('Working with tomography with different scaling for different layers.') self.tomography_params = Core_scaling_w_constrain() elif self.tomography_type == "mass_of_core_wo_constrain": - logging.debug('Working without any external constraints') + logging.debug('Working with tomography without any external constraints.') self.tomography_params = Core_scaling_wo_constrain() @@ -326,10 +375,8 @@ def calc_probs(self, nubar, e_array, rho_array, len_array, out): else: mix_matrix = self.osc_params.mix_matrix_complex - logging.debug('matter potential:\n%s' - % self.gen_mat_pot_matrix_complex) - logging.debug('decay matrix:\n%s' - % self.decay_matix) + logging.debug('matter potential:\n%s', self.gen_mat_pot_matrix_complex) + logging.debug('decay matrix:\n%s', self.decay_matix) propagate_array(self.osc_params.dm_matrix, # pylint: disable = unexpected-keyword-arg, no-value-for-parameter mix_matrix, @@ -430,24 +477,29 @@ def compute_function(self): # now we can proceed to calculate the generalised matter potential matrix std_mat_pot_matrix = np.zeros((3, 3), dtype=FTYPE) + 1.j * np.zeros((3, 3), dtype=FTYPE) - std_mat_pot_matrix[0, 0] += 1.0 + if not self.include_nlo: + logging.debug("Proceeding with *tree-level* standard CC matter potential....") + std_mat_pot_matrix[0, 0] += 1.0 + else: + logging.debug("Proceeding with *NLO* standard CC matter potential....") + std_mat_pot_matrix[0, 0] += 1.020 # add effective nsi coupling matrix if self.nsi_type is not None: - logging.debug('NSI matrix:\n%s' % self.nsi_params.eps_matrix) + logging.debug('NSI matrix:\n%s', self.nsi_params.eps_matrix) self.gen_mat_pot_matrix_complex = ( std_mat_pot_matrix + self.nsi_params.eps_matrix ) - logging.debug('Using generalised matter potential:\n%s' - % self.gen_mat_pot_matrix_complex) + logging.debug('Using generalised matter potential:\n%s', + self.gen_mat_pot_matrix_complex) else: self.gen_mat_pot_matrix_complex = std_mat_pot_matrix - logging.debug('Using standard matter potential:\n%s' - % self.gen_mat_pot_matrix_complex) + logging.debug('Using standard matter potential:\n%s', + self.gen_mat_pot_matrix_complex) if self.neutrino_decay: self.decay_matix = self.decay_params.decay_matrix - logging.debug('Decay matrix:\n%s' % self.decay_params.decay_matrix) + logging.debug('Decay matrix:\n%s', self.decay_params.decay_matrix) else : self.decay_matix = np.zeros((3, 3), dtype=FTYPE) + 1.j * np.zeros((3, 3), dtype=FTYPE) From 4cc521a1c5e4e347547a237f9f88d424881a82d6 Mon Sep 17 00:00:00 2001 From: T Ehrhardt Date: Wed, 10 Dec 2025 17:09:37 +0100 Subject: [PATCH 2/7] expand prob3 docstring (physics effects) + some cosmetic changes --- pisa/stages/osc/prob3.py | 153 +++++++++++++++++++++++++-------------- 1 file changed, 100 insertions(+), 53 deletions(-) diff --git a/pisa/stages/osc/prob3.py b/pisa/stages/osc/prob3.py index 7518f0d38..1cd8b3fd8 100644 --- a/pisa/stages/osc/prob3.py +++ b/pisa/stages/osc/prob3.py @@ -1,8 +1,6 @@ """ -PISA pi stage for the calculation of earth layers and osc. probabilities - -Maybe it would amke sense to split this up into a separate earth layer stage -and an osc. stage....todo +Service for the calculation of three-flavour oscillation probabilities, +allowing for various non-standard effects. """ @@ -35,42 +33,97 @@ class prob3(Stage): # pylint: disable=invalid-name """ - Extended Prob3-like oscillations class + Extended Prob3-like oscillations class. + + Expected container keys are: + "true_energy", + "true_coszen", + "nubar", + "flav", + "nu_flux", + "weights" Parameters ---------- include_nlo : bool (default: `False`) - Whether to include a +2.0% NLO correction to the SM CC - matter potential, as per - https://inspirehep.net/literature/2914951. + Whether to include a +2.0% NLO correction to the SM CC matter potential, + as per https://inspirehep.net/literature/2914951 (PRD111(2025)11). nsi_type : str or `None` (default: `None`) Choice of propagation/NC NSI parameterization. - If string, either 'standard' or 'vacuum-like'. + If string, either 'standard' or 'vacuum-like' + (see e.g. https://inspirehep.net/literature/1672932 (JHEP08(2018)180). + Parameters of the 'standard' parameterization are: + `eps_ee`, `eps_mumu`, `eps_tautau`, `eps_emu_magn`, `eps_emu_phase`, + `eps_etau_magn`, `eps_etau_phase`, `eps_mutau_magn`, `eps_mutau_phase`. + Parameters of the 'vacuum-like' one are: + `eps_scale`, `eps_prime`, `phi12`, `phi13`, `phi23`, `alpha1`, `alpha2`, + `deltansi`. reparam_mix_matrix : bool (default: `False`) - Whether to rephase the parameterization of the - leptonic mixing matrix from its PDG default by - diag(e^(i*delta_CP), 1, 1). This has no impact on - oscillation probabilities in the *absence* of NSI. - TODO: motivation + Whether to rephase the parameterization of the leptonic mixing matrix + from its PDG default by :math:`\mathrm{diag}(e^{i\delta_\mathrm{CP}}, 1, 1)`, + as motivated in https://inspirehep.net/literature/1672932 (JHEP08(2018)180). + In the *absence* of NSI, this has no observable impact on oscillation + probabilities, but results in the CPT transformation being realised + by the transformations :math:`\Delta m^2_{31} \\to -\Delta m^2_{32}, + \\theta_{12} \\to \pi/2 - \\theta_{12}, + \delta_\mathrm{CP} \\to \pi - \delta_\mathrm{CP}`, + which is accommodated by the parameters' usual ranges. These + transformations are then part of the generalized mass ordering + degeneracy in the presence of NSI (see same reference above). neutrino_decay : bool (default: `False`) - Whether to invoke neutrino decay with oscillations. - TODO: details + Whether to invoke neutrino decay with oscillations. The model + implemented assumes the (invisible) decay of the third mass eigenstate + into a singlet that is lighter than all three active mass eigenstates. + The decay is parameterized in a model-independent way via the parameter + :math:`\\alpha_3 = m_3/\\tau_3`, where :math:`m_3` is the mass and + :math:`\\tau_3` the lifetime in the rest frame. As a result, the vacuum + Hamiltonian in the flavor basis then reads + :math:`U \mathrm{diag}(0, \Delta m^2_{21}, \Delta m^2_{31} - i \\alpha_3) U^\dagger/(2E)`, + which deviates from the standard expression by the imaginary subtrahend in the + last entry. See e.g. https://inspirehep.net/literature/2870979 + (JHEP04(2025)105) for such an analysis in the literature. + The parameter :math:`\\alpha_3` is implemented by `decay_alpha3`. tomography_type : str or `None` (default: `None`) + Whether to allow for certain Earth matter density variations. If string, either 'mass_of_earth', 'mass_of_core_w_constrain', or 'mass_of_core_wo_constrain'. + In case of 'mass_of_earth', expects the single parameter + `density_scale`, which acts as an overall density scaling factor of + the assumed density profile. + In case of 'mass_of_core_w_constrain', expects the single parameter + `core_density_scale`, which acts as a density scaling factor for the + inner and outer core of the Earth in a 5-layer density model, + conserving the total mass and moment of inertia of the Earth by + accordingly rescaling the densities of the "inner" and "middle" (but + not the "outer") mantle. + In case of 'mass_of_core_wo_constrain', expects the three parameters + `core_density_scale`, `innermantle_density_scale`, + `middlemantle_density_scale`. In this case, the corresponding densities + in the 5-layer model are independently scaled, in contrast to the choice + above. + See e.g. https://inspirehep.net/literature/3072379 + (Eur.Phys.J.ST234(2025)16,5055-5064) for example analyses. lri_type : str or `None` (default: `None`) - Choice of model/parameterisation of long-range interactions. - If string, either 'emu-symmetry', 'etau-symmetry', - or 'mutau-symmetry'. + Choice of model/parameterization of long-range interactions (LRI). + If string, either 'emu-symmetry', 'etau-symmetry', or 'mutau-symmetry'. + Implemented is the single parameter `v_lri`. + In the case of :math:`e`-:math:`\mu` symmetry, it is added to the standard + :math:`ee` matter potential entry and subtracted from the :math:`\mu\mu` one. + In the case of :math:`e`-:math:`\\tau` symmetry, it is added to the standard + :math:`ee` matter potential entry and subtracted from the :math:`\\tau\\tau` one. + In the case of :math:`\mu`-:math:`\\tau` symmetry, it is added to the standard + :math:`\mu\mu` matter potential entry and subtracted from the :math:`\\tau\\tau` + one. See e.g. https://inspirehep.net/literature/2658147 (JHEP 08(2023)101) + for such an analysis in the literature. params - Expected params .. :: + expected params are .. :: detector_depth : float earth_model : PREM file path @@ -108,18 +161,8 @@ class prob3(Stage): # pylint: disable=invalid-name decay_alpha3 : quantity (energy^2) v_lri : quantity (eV) - Expected container keys are .. :: - - "true_energy" - "true_coszen" - "nubar" - "flav" - "nu_flux" - "weights" - - **kwargs + **std_kwargs Other kwargs are handled by Stage - ----- """ @@ -168,26 +211,25 @@ def __init__( nsi_type = nsi_type.strip().lower() if not nsi_type in choices: raise ValueError( - 'Chosen NSI type "%s" not available! Choose one of %s.' - % (nsi_type, choices) + f'Chosen NSI type "{nsi_type}" not available!' + f' Choose one of {choices}.' ) self.nsi_type = nsi_type """Type of NSI to assume.""" self.reparam_mix_matrix = reparam_mix_matrix """Use a PMNS mixing matrix parameterisation that differs from - the standard one by an overall phase matrix - diag(e^(i*delta_CP), 1, 1). This has no impact on - oscillation probabilities in the *absence* of NSI.""" + the standard one by an overall phase matrix + :math:`\mathrm{diag}(e^{i\delta_\mathrm{CP}}, 1, 1)`. This has no + impact on oscillation probabilities in the *absence* of NSI.""" self.neutrino_decay = neutrino_decay + """Invoke neutrino decay with neutrino oscillation.""" if neutrino_decay: self.decay_flag = 1 - else : + else: self.decay_flag = -1 - """Invoke neutrino decay with neutrino oscillation.""" - if self.nsi_type is None: nsi_params = () @@ -213,7 +255,7 @@ def __init__( 'eps_tautau' ) - if self.neutrino_decay : + if self.neutrino_decay: decay_params = ('decay_alpha3',) else: decay_params = () @@ -223,10 +265,11 @@ def __init__( lri_type = lri_type.strip().lower() if not lri_type in choices: raise ValueError( - 'Chosen LRI symmetry type "%s" not available! Choose one of %s.' - % (lri_type, choices) + f'Chosen LRI symmetry type "{lri_type}" not available!' + f' Choose one of {choices}.' ) self.lri_type = lri_type + """Type of LRI to assume.""" if self.lri_type is None: lri_params = () @@ -241,8 +284,8 @@ def __init__( choices = TOMOGRAPHY_TYPES if not tomography_type in choices: raise ValueError( - 'Chosen tomography type "%s" not available! Choose one of %s.' - % (tomography_type, choices) + f'Chosen tomography type "{tomography_type}" not available!' + f' Choose one of {choices}.' ) if tomography_type == 'mass_of_earth': tomography_params = ('density_scale',) @@ -254,6 +297,7 @@ def __init__( 'middlemantle_density_scale' ) self.tomography_type = tomography_type + """Type of Earth tomography to assume.""" expected_params = (expected_params + nsi_params + decay_params @@ -376,13 +420,13 @@ def calc_probs(self, nubar, e_array, rho_array, len_array, out): mix_matrix = self.osc_params.mix_matrix_complex logging.debug('matter potential:\n%s', self.gen_mat_pot_matrix_complex) - logging.debug('decay matrix:\n%s', self.decay_matix) + logging.debug('decay matrix:\n%s', self.decay_matrix) propagate_array(self.osc_params.dm_matrix, # pylint: disable = unexpected-keyword-arg, no-value-for-parameter mix_matrix, self.gen_mat_pot_matrix_complex, self.decay_flag, - self.decay_matix, + self.decay_matrix, self.lri_pot, nubar, e_array, @@ -406,7 +450,9 @@ def compute_function(self): YeM = self.params.YeM.value.m_as('dimensionless') if YeI != self.YeI or YeO != self.YeO or YeM != self.YeM: - self.YeI = YeI; self.YeO = YeO; self.YeM = YeM + self.YeI = YeI + self.YeO = YeO + self.YeM = YeM self.layers.setElecFrac(self.YeI, self.YeO, self.YeM) for container in self.data: self.layers.calcLayers(container['true_coszen']) @@ -417,8 +463,9 @@ def compute_function(self): # some safety checks on units # trying to avoid issue of angles with no dimension being assumed to be radians # here we enforce the user must speficy a valid angle unit - for angle_param in [self.params.theta12, self.params.theta13, self.params.theta23, self.params.deltacp] : - assert angle_param.value.units != ureg.dimensionless, "Param %s is dimensionless, but should have angle units [rad, degree]" % angle_param.name + for angle_param in [self.params.theta12, self.params.theta13, self.params.theta23, self.params.deltacp]: + if angle_param.value.units == ureg.dimensionless: + raise ValueError(f"{angle_param.name} is dimensionless, but needs units rad or deg!") # --- update mixing params --- self.osc_params.theta12 = self.params.theta12.value.m_as('rad') @@ -498,10 +545,10 @@ def compute_function(self): self.gen_mat_pot_matrix_complex) if self.neutrino_decay: - self.decay_matix = self.decay_params.decay_matrix + self.decay_matrix = self.decay_params.decay_matrix logging.debug('Decay matrix:\n%s', self.decay_params.decay_matrix) - else : - self.decay_matix = np.zeros((3, 3), dtype=FTYPE) + 1.j * np.zeros((3, 3), dtype=FTYPE) + else: + self.decay_matrix = np.zeros((3, 3), dtype=FTYPE) + 1.j * np.zeros((3, 3), dtype=FTYPE) self.lri_pot = np.zeros((3, 3), dtype=FTYPE) types_lri = ['emu-symmetry', 'etau-symmetry', 'etau-symmetry'] @@ -577,4 +624,4 @@ def init_test(**param_kwargs): Param(name='deltam31', value=3e-3*ureg.eV**2, **param_kwargs), Param(name='deltacp', value=180*ureg.degree, **param_kwargs), ]) - return prob3(params=param_set) + return prob3(include_nlo=True, params=param_set) From 4f42c4022ae30c3e1ae55d21efb1783577098a54 Mon Sep 17 00:00:00 2001 From: T Ehrhardt Date: Wed, 10 Dec 2025 22:20:23 +0100 Subject: [PATCH 3/7] clean up scaling_params module: boilerplate code, unnecessary imports, functions, variables; introduce globals for assumed 5-layer density profile; some cosmetics --- pisa/stages/osc/scaling_params.py | 162 +++++++++++------------------- 1 file changed, 60 insertions(+), 102 deletions(-) diff --git a/pisa/stages/osc/scaling_params.py b/pisa/stages/osc/scaling_params.py index 9cfac47e1..0f2b1c4af 100644 --- a/pisa/stages/osc/scaling_params.py +++ b/pisa/stages/osc/scaling_params.py @@ -3,130 +3,100 @@ Date : August 10, 2023 """ -from __future__ import absolute_import,division - import numpy as np -import os from pisa import FTYPE -import numba -if numba is None: - class jit(object): - """Decorator class to mimic Numba's `jit` when Numba is missing""" - def __init__(self, *args, **kwargs): - pass - def __call__(self, *args): - return args[0] -else: - jit = numba.jit - ftype = numba.typeof(FTYPE(1)) +__all__ = [ + 'Mass_scaling', 'Core_scaling_w_constrain', 'Core_scaling_wo_constrain', + 'FIVE_LAYER_RADII', 'FIVE_LAYER_RHOS' +] +FIVE_LAYER_RADII = np.array([0.0, 1221.50, 3480.00, 5701.00, 6151.0, 6371.00], dtype=FTYPE) +"""Radii (km) of five Earth layers to assume for the last two types of tomography.""" +FIVE_LAYER_RHOS = np.array([13.0, 13.0, 10.96, 5.03, 3.7, 2.5], dtype=FTYPE) +"""Matter densities (g/cm^3) of five Earth layers to assume for the last two types of tomography.""" -__all__ = ['Mass_scaling','Core_scaling_w_constrain','Core_scaling_wo_constrain'] class Mass_scaling(): """ - Uses a single scaling factor for all the layers. Scaling factor can be only positive. + Uses a single positive scaling factor for all the layers. """ def __init__(self): self._density_scale = 0. @property def density_scale(self): - return self._density_scale - + @density_scale.setter def density_scale(self, value): + assert value >= 0. self._density_scale = value - -class Core_scaling_w_constrain(object): + +class Core_scaling_w_constrain(): """ - Returns scaling factors for inner mantle and middle mantle by taking scaling factor of inner core and outer core as input. - Scaling factor of inner and outer core = core_density_scale (alpha) - Scaling factor of inner mantle = beta - Scaling factor of middle mantle = gamma - Outer mantle not scaled - This function solves the equations for two constraints: mass of earth and moment of inertia, by taking core_density_scale as an independent - parameter, and returns scaling factor factors for inner and middle mantle. - + Returns scaling factors for inner mantle and middle mantle by taking scaling + factor of inner core and outer core as input. + Scaling factor of inner and outer core = core_density_scale (alpha). + Scaling factor of inner mantle = beta. + Scaling factor of middle mantle = gamma. + Outer mantle not scaled. + This function solves the equations for two constraints: mass of earth and + moment of inertia, by taking core_density_scale as an independent + parameter, and returns scaling factors for inner and middle mantle. + """ def __init__(self): self._core_density_scale = 0. @property def core_density_scale(self): - return self._core_density_scale - + @core_density_scale.setter def core_density_scale(self, value): self._core_density_scale = value - def is_positive(self,lst): - for i in range(len(lst)): - if lst[i] < 0: - return False - return True - - def is_descending(self,lst): - for i in range(len(lst) - 1): - if lst[i] < lst[i + 1]: - return False - return True - @property def scaling_array(self): + radii_cm = FIVE_LAYER_RADII * 10**5 + rho = FIVE_LAYER_RHOS + + a1 = (4*np.pi/3)*(rho[1]*radii_cm[1]**3) + a2 = (8*np.pi/15)*(rho[1]*radii_cm[1]**5) + b1 = (4*np.pi/3)*(rho[2]*(radii_cm[2]**3 - radii_cm--[1]**3)) + b2 = (8*np.pi/15)*(rho[2]*(radii_cm[2]**5 - radii_cm[1]**5)) + c1 = (4*np.pi/3)*(rho[3]*(radii_cm[3]**3 - radii_cm[2]**3)) + c2 = (8*np.pi/15)*(rho[3]*(radii_cm[3]**5 - radii_cm[2]**5)) + d1 = (4*np.pi/3)*(rho[4]*(radii_cm[4]**3 - radii_cm[3]**3)) + d2 = (8*np.pi/15)*(rho[4]*(radii_cm[4]**5 - radii_cm[3]**5)) + e1 = (4*np.pi/3)*(rho[5]*(radii_cm[5]**3 - radii_cm[4]**3)) + e2 = (8*np.pi/15)*(rho[5]*(radii_cm[5]**5 - radii_cm[4]**5)) + + I = a2 + b2 + c2 + d2 + e2 + M = a1 + b1 + c1 + d1 + e1 - radius= np.array([0.0, 1221.50, 3480.00, 5701.00, 6151.0, 6371.00]) - R = radius * 10**5 - - rho = np.array([13.0, 13.0, 10.96, 5.03, 3.7, 2.5]) + alpha = self.core_density_scale + gamma = ((I*c1 - M*c2) - alpha*(c1*a2 - c2*a1) - alpha*(c1*b2 - b1*c2)-(c1*e2 - e1*c2))/(c1*d2 - d1*c2) + beta = (I - alpha*a2 - alpha*b2 - gamma*d2 - e2)/c2 - a1 = (4*np.pi/3)*(rho[1]* R[1]**3) - a2 = (8*np.pi/15)*(rho[1]* R[1]**5) - b1 = (4*np.pi/3)*(rho[2]* (R[2]**3 - R[1]**3)) - b2 = (8*np.pi/15)*(rho[2]* (R[2]**5 - R[1]**5)) - c1 = (4*np.pi/3)*(rho[3]* (R[3]**3 - R[2]**3)) - c2 = (8*np.pi/15)*(rho[3]* (R[3]**5 - R[2]**5)) - d1 = (4*np.pi/3)*(rho[4]* (R[4]**3 - R[3]**3)) - d2 = (8*np.pi/15)*(rho[4]* (R[4]**5 - R[3]**5)) - e1 = (4*np.pi/3)*(rho[5]* (R[5]**3 - R[4]**3)) - e2 = (8*np.pi/15)*(rho[5]* (R[5]**5 - R[4]**5)) + # density scaling factors need to be positive + assert (np.asarray([alpha, beta, gamma], dtype=FTYPE) >= 0).all() - I = a2 + b2 +c2 + d2 + e2 - M = a1 + b1 +c1 + d1 + e1 + tmp_array = np.ones(6, dtype=FTYPE) + tmp_array[1] = gamma + tmp_array[2] = beta + tmp_array[3] = alpha + tmp_array[4] = alpha + tmp_array[5] = alpha - alpha = self.core_density_scale - - - new_rho = np.zeros(6, dtype=FTYPE) - gamma = ((I*c1-M*c2)-alpha*(c1*a2 - c2*a1)- alpha*(c1*b2-b1*c2)-(c1*e2 - e1*c2))/(c1*d2-d1*c2) - beta = (I - alpha * a2 - alpha * b2 - gamma*d2 - e2)/(c2) - - - new_rho[0] = alpha * rho[0] - new_rho[1] = alpha * rho[1] - new_rho[2] = alpha * rho[2] - new_rho[3] = beta * rho[3] - new_rho[4] = gamma * rho[4] - new_rho[5] = rho[5] - - tmp_array = np.ones(6,dtype=FTYPE) - if self.is_positive(new_rho): # and self.is_descending(new_rho): ##turn this on if you want to put hydrostatic equilibrium condition - tmp_array[1] = gamma - tmp_array[2] = beta - tmp_array[3] = alpha - tmp_array[4] = alpha - tmp_array[5] = alpha - return tmp_array -class Core_scaling_wo_constrain(object): +class Core_scaling_wo_constrain(): """ - Takes scaling factors for core, inner mantle and outer mantle from pipeline and stores them in an array - + Stores independent scaling factors for core, inner mantle and outer mantle. + """ def __init__(self): self._core_density_scale = 0. @@ -135,55 +105,43 @@ def __init__(self): @property def core_density_scale(self): - return self._core_density_scale - + @core_density_scale.setter def core_density_scale(self, value): self._core_density_scale = value @property def innermantle_density_scale(self): - return self._innermantle_density_scale - + @innermantle_density_scale.setter def innermantle_density_scale(self, value): self._innermantle_density_scale = value @property def middlemantle_density_scale(self): - return self._middlemantle_density_scale - + @middlemantle_density_scale.setter def middlemantle_density_scale(self, value): self._middlemantle_density_scale = value @property def scaling_factor_array(self): - - tmp_array = np.ones(6,dtype=FTYPE) + tmp_array = np.ones(6, dtype=FTYPE) tmp_array[1] = self.middlemantle_density_scale tmp_array[2] = self.innermantle_density_scale tmp_array[3] = self.core_density_scale tmp_array[4] = self.core_density_scale tmp_array[5] = self.core_density_scale - + return tmp_array def test_scaling_params(): + #TODO pass if __name__=='__main__': - from pisa import TARGET - from pisa.utils.log import set_verbosity, logging - assert TARGET == 'cpu', "Cannot test functions on GPU, set PISA_TARGET to 'cpu'" - set_verbosity(1) test_scaling_params() - - - - - From 41bf5b69ae02449f2a30bdc06ef97e0b5efcbd60 Mon Sep 17 00:00:00 2001 From: T Ehrhardt Date: Thu, 11 Dec 2025 02:00:00 +0100 Subject: [PATCH 4/7] External vs. internal Earth model consistency check for tomography types where the five layers are required. --- pisa/stages/osc/prob3.py | 66 ++++++++++++++++++++++++---------------- 1 file changed, 39 insertions(+), 27 deletions(-) diff --git a/pisa/stages/osc/prob3.py b/pisa/stages/osc/prob3.py index 1cd8b3fd8..f61216204 100644 --- a/pisa/stages/osc/prob3.py +++ b/pisa/stages/osc/prob3.py @@ -16,7 +16,10 @@ from pisa.stages.osc.osc_params import OscParams from pisa.stages.osc.decay_params import DecayParams from pisa.stages.osc.lri_params import LRIParams -from pisa.stages.osc.scaling_params import Mass_scaling, Core_scaling_w_constrain, Core_scaling_wo_constrain +from pisa.stages.osc.scaling_params import ( + Mass_scaling, Core_scaling_w_constrain, Core_scaling_wo_constrain, + FIVE_LAYER_RADII, FIVE_LAYER_RHOS +) from pisa.stages.osc.layers import Layers from pisa.stages.osc.prob3numba.numba_osc_hostfuncs import propagate_array, fill_probs from pisa.utils.resources import find_resource @@ -32,7 +35,7 @@ class prob3(Stage): # pylint: disable=invalid-name - """ + r""" Extended Prob3-like oscillations class. Expected container keys are: @@ -67,9 +70,9 @@ class prob3(Stage): # pylint: disable=invalid-name as motivated in https://inspirehep.net/literature/1672932 (JHEP08(2018)180). In the *absence* of NSI, this has no observable impact on oscillation probabilities, but results in the CPT transformation being realised - by the transformations :math:`\Delta m^2_{31} \\to -\Delta m^2_{32}, - \\theta_{12} \\to \pi/2 - \\theta_{12}, - \delta_\mathrm{CP} \\to \pi - \delta_\mathrm{CP}`, + by the transformations :math:`\Delta m^2_{31} \to -\Delta m^2_{32}, + \theta_{12} \to \pi/2 - \theta_{12}, + \delta_\mathrm{CP} \to \pi - \delta_\mathrm{CP}`, which is accommodated by the parameters' usual ranges. These transformations are then part of the generalized mass ordering degeneracy in the presence of NSI (see same reference above). @@ -79,14 +82,14 @@ class prob3(Stage): # pylint: disable=invalid-name implemented assumes the (invisible) decay of the third mass eigenstate into a singlet that is lighter than all three active mass eigenstates. The decay is parameterized in a model-independent way via the parameter - :math:`\\alpha_3 = m_3/\\tau_3`, where :math:`m_3` is the mass and - :math:`\\tau_3` the lifetime in the rest frame. As a result, the vacuum + :math:`\alpha_3 = m_3/\tau_3`, where :math:`m_3` is the mass and + :math:`\tau_3` the lifetime in the rest frame. As a result, the vacuum Hamiltonian in the flavor basis then reads - :math:`U \mathrm{diag}(0, \Delta m^2_{21}, \Delta m^2_{31} - i \\alpha_3) U^\dagger/(2E)`, + :math:`U \mathrm{diag}(0, \Delta m^2_{21}, \Delta m^2_{31} - i \alpha_3) U^\dagger/(2E)`, which deviates from the standard expression by the imaginary subtrahend in the last entry. See e.g. https://inspirehep.net/literature/2870979 (JHEP04(2025)105) for such an analysis in the literature. - The parameter :math:`\\alpha_3` is implemented by `decay_alpha3`. + The parameter :math:`\alpha_3` is implemented by `decay_alpha3`. tomography_type : str or `None` (default: `None`) Whether to allow for certain Earth matter density variations. @@ -97,7 +100,7 @@ class prob3(Stage): # pylint: disable=invalid-name the assumed density profile. In case of 'mass_of_core_w_constrain', expects the single parameter `core_density_scale`, which acts as a density scaling factor for the - inner and outer core of the Earth in a 5-layer density model, + inner and outer core of the Earth in a pre-determined 5-layer density model, conserving the total mass and moment of inertia of the Earth by accordingly rescaling the densities of the "inner" and "middle" (but not the "outer") mantle. @@ -115,10 +118,10 @@ class prob3(Stage): # pylint: disable=invalid-name Implemented is the single parameter `v_lri`. In the case of :math:`e`-:math:`\mu` symmetry, it is added to the standard :math:`ee` matter potential entry and subtracted from the :math:`\mu\mu` one. - In the case of :math:`e`-:math:`\\tau` symmetry, it is added to the standard - :math:`ee` matter potential entry and subtracted from the :math:`\\tau\\tau` one. - In the case of :math:`\mu`-:math:`\\tau` symmetry, it is added to the standard - :math:`\mu\mu` matter potential entry and subtracted from the :math:`\\tau\\tau` + In the case of :math:`e`-:math:`\tau` symmetry, it is added to the standard + :math:`ee` matter potential entry and subtracted from the :math:`\tau\tau` one. + In the case of :math:`\mu`-:math:`\tau` symmetry, it is added to the standard + :math:`\mu\mu` matter potential entry and subtracted from the :math:`\tau\tau` one. See e.g. https://inspirehep.net/literature/2658147 (JHEP 08(2023)101) for such an analysis in the literature. @@ -218,7 +221,7 @@ def __init__( """Type of NSI to assume.""" self.reparam_mix_matrix = reparam_mix_matrix - """Use a PMNS mixing matrix parameterisation that differs from + r"""Use a PMNS mixing matrix parameterisation that differs from the standard one by an overall phase matrix :math:`\mathrm{diag}(e^{i\delta_\mathrm{CP}}, 1, 1)`. This has no impact on oscillation probabilities in the *absence* of NSI.""" @@ -356,18 +359,6 @@ def setup_function(self): logging.debug('Working with LRI') self.lri_params = LRIParams() - - if self.tomography_type == "mass_of_earth": - logging.debug('Working with tomography with a single density scaling factor.') - self.tomography_params = Mass_scaling() - elif self.tomography_type == "mass_of_core_w_constrain": - logging.debug('Working with tomography with different scaling for different layers.') - self.tomography_params = Core_scaling_w_constrain() - elif self.tomography_type == "mass_of_core_wo_constrain": - logging.debug('Working with tomography without any external constraints.') - self.tomography_params = Core_scaling_wo_constrain() - - # setup the layers #if self.params.earth_model.value is not None: earth_model = find_resource(self.params.earth_model.value) @@ -379,6 +370,27 @@ def setup_function(self): self.layers = Layers(earth_model, detector_depth, prop_height) self.layers.setElecFrac(self.YeI, self.YeO, self.YeM) + if self.tomography_type == "mass_of_earth": + logging.debug('Working with tomography with a single density scaling factor.') + self.tomography_params = Mass_scaling() + else: + # not elegant but safe: external Earth model must correspond to internal hard-coded one + assert self.layers.using_earth_model + radii_equal = np.allclose(np.add(self.layers.radii[::-1][:-1], 1), np.add(FIVE_LAYER_RADII, 1)) + rhos_equal = np.allclose(np.add(self.layers.rhos_unweighted[::-1][:-1], 1), np.add(FIVE_LAYER_RHOS, 1)) + if not radii_equal or not rhos_equal: + raise ValueError( + "The Earth model provided needs to have the same layer radii" + " and densities as the one hard-coded for the chosen type of" + f" tomography internally, namely radii {FIVE_LAYER_RADII} and" + f" densities {FIVE_LAYER_RHOS}." + ) + if self.tomography_type == "mass_of_core_w_constrain": + logging.debug('Working with tomography with different scaling for different layers.') + self.tomography_params = Core_scaling_w_constrain() + elif self.tomography_type == "mass_of_core_wo_constrain": + logging.debug('Working with tomography without any external constraints.') + self.tomography_params = Core_scaling_wo_constrain() # --- calculate the layers --- if self.is_map: From 1688f2924048ab1da37303169323093e09bd9b83 Mon Sep 17 00:00:00 2001 From: T Ehrhardt Date: Thu, 11 Dec 2025 02:11:42 +0100 Subject: [PATCH 5/7] had accidentally applied checks even when tomography_type = None --- pisa/stages/osc/prob3.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/pisa/stages/osc/prob3.py b/pisa/stages/osc/prob3.py index f61216204..56f497b6d 100644 --- a/pisa/stages/osc/prob3.py +++ b/pisa/stages/osc/prob3.py @@ -373,11 +373,14 @@ def setup_function(self): if self.tomography_type == "mass_of_earth": logging.debug('Working with tomography with a single density scaling factor.') self.tomography_params = Mass_scaling() - else: + elif self.tomography_type is not None: # not elegant but safe: external Earth model must correspond to internal hard-coded one assert self.layers.using_earth_model - radii_equal = np.allclose(np.add(self.layers.radii[::-1][:-1], 1), np.add(FIVE_LAYER_RADII, 1)) - rhos_equal = np.allclose(np.add(self.layers.rhos_unweighted[::-1][:-1], 1), np.add(FIVE_LAYER_RHOS, 1)) + radii_ext = self.layers.radii[::-1][:-1] + rhos_ext = self.layers.rhos_unweighted[::-1][:-1] + assert len(radii_ext) == len(FIVE_LAYER_RADII) and len(rhos_ext) == len(FIVE_LAYER_RHOS) + radii_equal = np.allclose(np.add(radii_ext, 1), np.add(FIVE_LAYER_RADII, 1)) + rhos_equal = np.allclose(np.add(rhos_ext, 1), np.add(FIVE_LAYER_RHOS, 1)) if not radii_equal or not rhos_equal: raise ValueError( "The Earth model provided needs to have the same layer radii" From 841b940eefc2bb2a24e3c4603a08e9aa17d20f15 Mon Sep 17 00:00:00 2001 From: T Ehrhardt Date: Thu, 11 Dec 2025 02:25:49 +0100 Subject: [PATCH 6/7] remove spurious dashes --- pisa/stages/osc/scaling_params.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pisa/stages/osc/scaling_params.py b/pisa/stages/osc/scaling_params.py index 0f2b1c4af..2467c00c2 100644 --- a/pisa/stages/osc/scaling_params.py +++ b/pisa/stages/osc/scaling_params.py @@ -65,7 +65,7 @@ def scaling_array(self): a1 = (4*np.pi/3)*(rho[1]*radii_cm[1]**3) a2 = (8*np.pi/15)*(rho[1]*radii_cm[1]**5) - b1 = (4*np.pi/3)*(rho[2]*(radii_cm[2]**3 - radii_cm--[1]**3)) + b1 = (4*np.pi/3)*(rho[2]*(radii_cm[2]**3 - radii_cm[1]**3)) b2 = (8*np.pi/15)*(rho[2]*(radii_cm[2]**5 - radii_cm[1]**5)) c1 = (4*np.pi/3)*(rho[3]*(radii_cm[3]**3 - radii_cm[2]**3)) c2 = (8*np.pi/15)*(rho[3]*(radii_cm[3]**5 - radii_cm[2]**5)) @@ -78,7 +78,7 @@ def scaling_array(self): M = a1 + b1 + c1 + d1 + e1 alpha = self.core_density_scale - gamma = ((I*c1 - M*c2) - alpha*(c1*a2 - c2*a1) - alpha*(c1*b2 - b1*c2)-(c1*e2 - e1*c2))/(c1*d2 - d1*c2) + gamma = ((I*c1 - M*c2) - alpha*(c1*a2 - c2*a1) - alpha*(c1*b2 - b1*c2) - (c1*e2 - e1*c2))/(c1*d2 - d1*c2) beta = (I - alpha*a2 - alpha*b2 - gamma*d2 - e2)/c2 # density scaling factors need to be positive From f2be139cad07fadd83c93af83c06298d4ed4a7cb Mon Sep 17 00:00:00 2001 From: T Ehrhardt Date: Thu, 11 Dec 2025 18:40:23 +0100 Subject: [PATCH 7/7] improve tomography exception handling slightly --- pisa/stages/osc/prob3.py | 22 +++++++++++++--------- pisa/stages/osc/scaling_params.py | 8 +++++++- 2 files changed, 20 insertions(+), 10 deletions(-) diff --git a/pisa/stages/osc/prob3.py b/pisa/stages/osc/prob3.py index 56f497b6d..87706f4ce 100644 --- a/pisa/stages/osc/prob3.py +++ b/pisa/stages/osc/prob3.py @@ -18,7 +18,7 @@ from pisa.stages.osc.lri_params import LRIParams from pisa.stages.osc.scaling_params import ( Mass_scaling, Core_scaling_w_constrain, Core_scaling_wo_constrain, - FIVE_LAYER_RADII, FIVE_LAYER_RHOS + FIVE_LAYER_RADII, FIVE_LAYER_RHOS, TOMOGRAPHY_ERROR_MSG ) from pisa.stages.osc.layers import Layers from pisa.stages.osc.prob3numba.numba_osc_hostfuncs import propagate_array, fill_probs @@ -371,23 +371,27 @@ def setup_function(self): self.layers.setElecFrac(self.YeI, self.YeO, self.YeM) if self.tomography_type == "mass_of_earth": + if not self.layers.using_earth_model: + # already fail here instead of in compute upon attempt to scale + # (to be on the safe side; not sure if this case can even happen) + raise ValueError( + "You need to provide some Earth model, whose densities can" + " be rescaled!" + ) logging.debug('Working with tomography with a single density scaling factor.') self.tomography_params = Mass_scaling() elif self.tomography_type is not None: # not elegant but safe: external Earth model must correspond to internal hard-coded one - assert self.layers.using_earth_model + if not self.layers.using_earth_model: + raise ValueError(TOMOGRAPHY_ERROR_MSG) radii_ext = self.layers.radii[::-1][:-1] rhos_ext = self.layers.rhos_unweighted[::-1][:-1] - assert len(radii_ext) == len(FIVE_LAYER_RADII) and len(rhos_ext) == len(FIVE_LAYER_RHOS) + if not (len(radii_ext) == len(FIVE_LAYER_RADII) and len(rhos_ext) == len(FIVE_LAYER_RHOS)): + raise ValueError(TOMOGRAPHY_ERROR_MSG) radii_equal = np.allclose(np.add(radii_ext, 1), np.add(FIVE_LAYER_RADII, 1)) rhos_equal = np.allclose(np.add(rhos_ext, 1), np.add(FIVE_LAYER_RHOS, 1)) if not radii_equal or not rhos_equal: - raise ValueError( - "The Earth model provided needs to have the same layer radii" - " and densities as the one hard-coded for the chosen type of" - f" tomography internally, namely radii {FIVE_LAYER_RADII} and" - f" densities {FIVE_LAYER_RHOS}." - ) + raise ValueError(TOMOGRAPHY_ERROR_MSG) if self.tomography_type == "mass_of_core_w_constrain": logging.debug('Working with tomography with different scaling for different layers.') self.tomography_params = Core_scaling_w_constrain() diff --git a/pisa/stages/osc/scaling_params.py b/pisa/stages/osc/scaling_params.py index 2467c00c2..5ecb62b85 100644 --- a/pisa/stages/osc/scaling_params.py +++ b/pisa/stages/osc/scaling_params.py @@ -8,14 +8,20 @@ __all__ = [ 'Mass_scaling', 'Core_scaling_w_constrain', 'Core_scaling_wo_constrain', - 'FIVE_LAYER_RADII', 'FIVE_LAYER_RHOS' + 'FIVE_LAYER_RADII', 'FIVE_LAYER_RHOS', 'TOMOGRAPHY_ERROR_MSG' ] FIVE_LAYER_RADII = np.array([0.0, 1221.50, 3480.00, 5701.00, 6151.0, 6371.00], dtype=FTYPE) """Radii (km) of five Earth layers to assume for the last two types of tomography.""" + FIVE_LAYER_RHOS = np.array([13.0, 13.0, 10.96, 5.03, 3.7, 2.5], dtype=FTYPE) """Matter densities (g/cm^3) of five Earth layers to assume for the last two types of tomography.""" +TOMOGRAPHY_ERROR_MSG = ("You need to provide the appropriate 5-layer Earth model," +f" which has the same layer radii ({FIVE_LAYER_RADII}) and densities ({FIVE_LAYER_RHOS})" +" as the one hard-coded for the chosen type of tomography internally.") +"""An error message simultaneously explaining the state of the code.""" + class Mass_scaling(): """