Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Flavor Transformations plus Earth matter #258

Open
wants to merge 172 commits into
base: main
Choose a base branch
from

Conversation

jpkneller
Copy link
Contributor

@jpkneller jpkneller commented Jun 16, 2023

A rewrite of the flavor transformation classes.

  1. Includes a new method called get_probabilities which returns a matrix of transition probabilities. This method is meant t replace the 8 prob_xy methods in the previous versions of the classes
  2. To avoid duplicating code, there are now two new intermediate classes called ThreeFlavorTransformation and FourFlavorTransformation between the abstract base class FlavorTransformation and the various flavor transformation classes.
  3. The ThreeFlavorTransformation class has two methods: one called vacuum and the other Earth_matter. These methods return the elements of the D matrix previously calculated in the flavor transformation classes. The Earth_matter method does an earth-matter effect calculation using the Sqa3Earth module. Compiling that module is described elsewhere. The FourFlavorTransformation class only has the vacuum method.
  4. The flavor transformation classes now have an optional argument to the constructor called AltAz. This is an astropy object which gives the altitude and azimuth of the supernova relative to the observer / detector. If not None, the earth matter effect will be done (if the altitude is less than zero)
  5. The previous version of the flavor transformation classes had two optional arguments to the constructor called mixing_angles and mass hierarchy. These have been replaced with the optional argument called mix_params which must be an instance of either the MixingParamters3Flavor or MixingParamters4Flavor classes in neutrino.py. If None, the default is to the use the normal mass hierarchy as before.

jpkneller added 30 commits June 16, 2023 09:42
A rewrite of the flavor transformation classes. 
1) Includes a new method called get_probabilities which returns a matrix of transition probabilities. This method is meant t replace the 8 prob_xy methods in the previous versions of the classes
2) To avoid duplicating code, there are now two new intermediate classes called ThreeFlavorTransformation and FourFlavorTransformation between the abstract base class FlavorTransformation and the various flavor transformation classes.
3) The ThreeFlavorTransformation class has two methods: one called vacuum and the other Earth_matter. These methods return the elements of the D matrix previously calculated in the flavor transformation classes. The Earth_matter method does an earth-matter effect calculation using the Sqa3Earth module. Compiling that module is described elsewhere. The FourFlavorTransformation class only has the vacuum method. 
4) The flavor transformation classes now have an optional argument to the constructor called AltAz. This is an astropy object which gives the altitude and azimuth of the supernova relative to the observer / detector. If not None, the earth matter effect will be done (if the altitude is less than zero)
5) The previous version of the flavor transformation classes had two optional arguments to the constructor called mixing_angles and mass hierarchy. These have been replaced with the optional argument called mix_params which must be an instance of either the MixingParamters3Flavor or MixingParamters4Flavor classes in neutrino.py. If None, the default is to the use the normal mass hierarchy as before.
…on class inputs

The transformation_type string argument to the generate_* functions has been replaced with a flavor_transformation argument which can be either a string or a flavor transformation class. If it is a string the get_transformation function is called to identify the class.
Added __str__ methods
Now uses the get_probabilities method of the flavor transformation classes rather than the prob_xy methods.
…newpy/SQA/src/mstl/physics/units and constants.cpp
…newpy/SQA/src/mstl/math2/spline/discontinuous.cpp
…py/SQA/src/mstl/math2/analysis/algorithm3.cpp
… python/snewpy/SQA/src/mstl/math2/algebra/column and row vectors.cpp
@jpkneller
Copy link
Contributor Author

The Earth-matter effect module is now a separate repo and is not required for SNEWPY to work.

@jpkneller jpkneller added the enhancement New feature or request label Nov 17, 2023
@jpkneller jpkneller added this to the v2.0 milestone Nov 17, 2023
Copy link
Contributor

@Sheshuk Sheshuk left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi Jim, thank you for updating this!
Having EMEWS decoupled makes it much easier for users to explore the SNEWPY without having to run additional steps (compilation, installation etc), in case they don't need this functionality.

Also it's great to finally have the transformation matrices in place!

I have two suggestions:

  1. As you can see, the tests are failing, because they try to use mh parameter in the constructors of the FlavorTransformations. So for backward compatibility we should still use mh input, if it is provided (and add a note in the code to clean it up for SNEWPY v2.0)
  2. Move the repeating code, which fills the transformation matrix from p_ee, pbar_ee... to a separate function (I added an example in comments) and use it. Having less repeating code means less space for possible copy-paste errors.

Comment on lines +419 to +431
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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These lines (filling the matrix) are repeated in each transformation class, which makes a lot copy-paste.

I suggest implementing matrix creation as a separate utility function, something like this:

def _gen_matrix(p_ee, p_ex=None, pbar_ee=None, pbar_ex=None):
    "Generate transformation matrix, given the probabilities"
    p_ee = np.array(p_ee)
    p_ex = np.array(p_ex) if (p_ex is not None) else 1-p_ee
    pbar_ee = np.array(pbar_ee) if(pbar_ee is not None) else p_ee
    pbar_ex = np.array(pbar_ex) if(pbar_ex is not None) else 1-pbar_ee
    
    p = np.zeros(shape=[4,4,*p_ee.shape])
    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

    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

and put it in the global scope in this file.

So this (and other places like this) will become just

Suggested change
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
return _gen_matrix(p_ee=p_ee, pbar_ee=pbar_ee)

Comment on lines +460 to +472
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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as above (use _gen_matrix function):

Suggested change
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
return _gen_matrix(p_ee=p_ee, pbar_ee=pbar_ee)

Comment on lines +518 to +530
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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as above (use _gen_matrix function):

Suggested change
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
return _gen_matrix(p_ee=p_ee, pbar_ee=pbar_ee)

Comment on lines +572 to +584
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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as above (use _gen_matrix function):

Suggested change
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
return _gen_matrix(p_ee=p_ee, pbar_ee=pbar_ee)

Comment on lines +629 to +641
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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as above (use _gen_matrix function):

Suggested change
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
return _gen_matrix(p_ee=p_ee, pbar_ee=pbar_ee)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These "suggestions" are looking weird, looks like it grabs the text from both panels ("main" and this branch). I hope it's still clear what I meant

Comment on lines +670 to +682
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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as above (use _gen_matrix function):

        return _gen_matrix(p_ee=p_ee, pbar_ee=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_ex/2

return p
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as above (use _gen_matrix function):

        return _gen_matrix(p_ee, p_ex, pbar_ee, pbar_ex)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants