-
Notifications
You must be signed in to change notification settings - Fork 22
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
base: main
Are you sure you want to change the base?
Conversation
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.
…/adiabatic_basis.cpp
…newpy/SQA/src/mstl/physics/units and constants.cpp
…stl/math2/number.h
…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
…to python/snewpy/EMEWS/src/mstl/math2/algebra/column and row vectors.cpp
…ewpy/EMEWS/src/mstl/math2/analysis/algorithm3.cpp
…/snewpy/EMEWS/src/mstl/math2/spline/discontinuous.cpp
…/snewpy/EMEWS/src/mstl/physics/units and constants.cpp
Moved EMEWS to separate repository
EMEWS is now optional
The Earth-matter effect module is now a separate repo and is not required for SNEWPY to work. |
Resolving conflicts in flux.py and snowglobes.py
There was a problem hiding this 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:
- 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 usemh
input, if it is provided (and add a note in the code to clean it up for SNEWPY v2.0) - 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.
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 |
There was a problem hiding this comment.
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
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) |
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 |
There was a problem hiding this comment.
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):
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) |
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 |
There was a problem hiding this comment.
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):
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) |
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 |
There was a problem hiding this comment.
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):
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) |
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 |
There was a problem hiding this comment.
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):
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) |
There was a problem hiding this comment.
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
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 |
There was a problem hiding this comment.
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 |
There was a problem hiding this comment.
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)
A rewrite of the flavor transformation classes.