aboutsummaryrefslogtreecommitdiffstats
path: root/utils
diff options
context:
space:
mode:
Diffstat (limited to 'utils')
-rw-r--r--utils/fr.py101
-rw-r--r--utils/gf.py1
-rw-r--r--utils/misc.py5
-rw-r--r--utils/multinest.py1
-rw-r--r--utils/plot.py2
5 files changed, 77 insertions, 33 deletions
diff --git a/utils/fr.py b/utils/fr.py
index f98c730..51864a3 100644
--- a/utils/fr.py
+++ b/utils/fr.py
@@ -14,14 +14,61 @@ from functools import partial
import numpy as np
from utils.enums import EnergyDependance
-from utils.misc import DTYPE, CDTYPE, PI
from utils.misc import enum_parse, parse_bool
+import mpmath as mp
+mp.mp.dps = 100 # Computation precision
+
+# DTYPE = np.float128
+# CDTYPE = np.complex256
+# PI = np.arccos(DTYPE(-1))
+# SQRT = np.sqrt
+# COS = np.cos
+# SIN = np.sin
+# ACOS = np.arccos
+# ASIN = np.arcsin
+# EXP = np.exp
+
+DTYPE = mp.mpf
+CDTYPE = mp.mpc
+PI = mp.pi
+SQRT = mp.sqrt
+COS = mp.cos
+SIN = mp.sin
+ACOS = mp.acos
+ASIN = mp.asin
+EXP = mp.exp
MASS_EIGENVALUES = [7.40E-23, 2.515E-21]
"""SM mass eigenvalues"""
+def determinant(x):
+ """Calculate the determininant of a 3x3 matrix.
+
+ Parameters
+ ----------
+ x : ndarray, shape = (3, 3)
+
+ Returns
+ ----------
+ float determinant
+
+ Examples
+ ----------
+ >>> print determinant(
+ >>> [[-1.65238188-0.59549718j, 0.27486548-0.18437467j, -1.35524534-0.38542072j],
+ >>> [-1.07480906+0.29630449j, -0.47808456-0.80316821j, -0.88609356-1.50737308j],
+ >>> [-0.14924144-0.99230446j, 0.49504234+0.63639805j, 2.29258915-0.36537507j]]
+ >>> )
+ (2.7797571563274688+3.0841795325804848j)
+
+ """
+ return (x[0][0] * (x[1][1] * x[2][2] - x[2][1] * x[1][2])
+ -x[1][0] * (x[0][1] * x[2][2] - x[2][1] * x[0][2])
+ +x[2][0] * (x[0][1] * x[1][2] - x[1][1] * x[0][2]))
+
+
def angles_to_fr(src_angles):
"""Convert angular projection of the source flavour ratio back into the
flavour ratio.
@@ -41,13 +88,13 @@ def angles_to_fr(src_angles):
(0.38340579025361626, 0.16431676725154978, 0.45227744249483393)
"""
- sphi4, c2psi = src_angles
+ sphi4, c2psi = map(DTYPE, src_angles)
- psi = (0.5)*np.arccos(c2psi, dtype=DTYPE)
+ psi = (0.5)*ACOS(c2psi)
- sphi2 = np.sqrt(sphi4, dtype=DTYPE)
+ sphi2 = SQRT(sphi4)
cphi2 = 1. - sphi2
- spsi2 = np.sin(psi, dtype=DTYPE)**2
+ spsi2 = SIN(psi)**2
cspi2 = 1. - spsi2
x = sphi2*cspi2
@@ -78,27 +125,27 @@ def angles_to_u(bsm_angles):
[ 0.28614067-0.42427084j, -0.64749908-0.21213542j, 0.52331757+0.j ]])
"""
- s12_2, c13_4, s23_2, dcp = bsm_angles
+ s12_2, c13_4, s23_2, dcp = map(DTYPE, bsm_angles)
dcp = CDTYPE(dcp)
c12_2 = 1. - s12_2
- c13_2 = np.sqrt(c13_4, dtype=DTYPE)
+ c13_2 = SQRT(c13_4)
s13_2 = 1. - c13_2
c23_2 = 1. - s23_2
- t12 = np.arcsin(np.sqrt(s12_2, dtype=DTYPE))
- t13 = np.arccos(np.sqrt(c13_2, dtype=DTYPE))
- t23 = np.arcsin(np.sqrt(s23_2, dtype=DTYPE))
+ t12 = ASIN(SQRT(s12_2))
+ t13 = ACOS(SQRT(c13_2))
+ t23 = ASIN(SQRT(s23_2))
- c12 = np.cos(t12)
- s12 = np.sin(t12)
- c13 = np.cos(t13)
- s13 = np.sin(t13)
- c23 = np.cos(t23)
- s23 = np.sin(t23)
+ c12 = COS(t12)
+ s12 = SIN(t12)
+ c13 = COS(t13)
+ s13 = SIN(t13)
+ c23 = COS(t23)
+ s23 = SIN(t23)
p1 = np.array([[1 , 0 , 0] , [0 , c23 , s23] , [0 , -s23 , c23]] , dtype=CDTYPE)
- p2 = np.array([[c13 , 0 , s13*np.exp(-1j*dcp)] , [0 , 1 , 0] , [-s13*np.exp(1j*dcp) , 0 , c13]] , dtype=CDTYPE)
+ p2 = np.array([[c13 , 0 , s13*EXP(-1j*dcp)] , [0 , 1 , 0] , [-s13*EXP(1j*dcp) , 0 , c13]] , dtype=CDTYPE)
p3 = np.array([[c12 , s12 , 0] , [-s12 , c12 , 0] , [0 , 0 , 1]] , dtype=CDTYPE)
u = np.dot(np.dot(p1, p2), p3)
@@ -141,15 +188,15 @@ def cardano_eqn(ham):
a = -np.trace(ham)
b = DTYPE(1)/2 * ((np.trace(ham))**DTYPE(2) - np.trace(np.dot(ham, ham)))
- c = -np.linalg.det(ham)
+ c = -determinant(ham)
Q = (DTYPE(1)/9) * (a**DTYPE(2) - DTYPE(3)*b)
R = (DTYPE(1)/54) * (DTYPE(2)*a**DTYPE(3) - DTYPE(9)*a*b + DTYPE(27)*c)
- theta = np.arccos(R / np.sqrt(Q**DTYPE(3)))
+ theta = ACOS(R / SQRT(Q**DTYPE(3)))
- E1 = -DTYPE(2) * np.sqrt(Q) * np.cos(theta/DTYPE(3)) - (DTYPE(1)/3)*a
- E2 = -DTYPE(2) * np.sqrt(Q) * np.cos((theta - DTYPE(2)*PI)/DTYPE(3)) - (DTYPE(1)/3)*a
- E3 = -DTYPE(2) * np.sqrt(Q) * np.cos((theta + DTYPE(2)*PI)/DTYPE(3)) - (DTYPE(1)/3)*a
+ E1 = -DTYPE(2) * SQRT(Q) * COS(theta/DTYPE(3)) - (DTYPE(1)/3)*a
+ E2 = -DTYPE(2) * SQRT(Q) * COS((theta - DTYPE(2)*PI)/DTYPE(3)) - (DTYPE(1)/3)*a
+ E3 = -DTYPE(2) * SQRT(Q) * COS((theta + DTYPE(2)*PI)/DTYPE(3)) - (DTYPE(1)/3)*a
A1 = ham[1][2] * (ham[0][0] - E1) - ham[1][0]*ham[0][2]
A2 = ham[1][2] * (ham[0][0] - E2) - ham[1][0]*ham[0][2]
@@ -163,9 +210,9 @@ def cardano_eqn(ham):
C2 = ham[1][0] * (ham[2][2] - E2) - ham[1][2]*ham[2][0]
C3 = ham[1][0] * (ham[2][2] - E3) - ham[1][2]*ham[2][0]
- N1 = np.sqrt(np.abs(A1*B1, dtype=DTYPE)**2 + np.abs(A1*C1, dtype=DTYPE)**2 + np.abs(B1*C1, dtype=DTYPE)**2)
- N2 = np.sqrt(np.abs(A2*B2, dtype=DTYPE)**2 + np.abs(A2*C2, dtype=DTYPE)**2 + np.abs(B2*C2, dtype=DTYPE)**2)
- N3 = np.sqrt(np.abs(A3*B3, dtype=DTYPE)**2 + np.abs(A3*C3, dtype=DTYPE)**2 + np.abs(B3*C3, dtype=DTYPE)**2)
+ N1 = SQRT(np.abs(A1*B1)**2 + np.abs(A1*C1)**2 + np.abs(B1*C1)**2)
+ N2 = SQRT(np.abs(A2*B2)**2 + np.abs(A2*C2)**2 + np.abs(B2*C2)**2)
+ N3 = SQRT(np.abs(A3*B3)**2 + np.abs(A3*C3)**2 + np.abs(B3*C3)**2)
mm = np.array([
[np.conjugate(B1)*C1 / N1, np.conjugate(B2)*C2 / N2, np.conjugate(B3)*C3 / N3],
@@ -383,7 +430,7 @@ def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES,
mass_matrix = np.array(
[[0, 0, 0], [0, mass_eigenvalues[0], 0], [0, 0, mass_eigenvalues[1]]]
)
- sm_ham = (1./(2*energy))*np.matmul(sm_u, np.matmul(mass_matrix, sm_u.conj().T))
+ sm_ham = (1./(2*energy))*np.dot(sm_u, np.dot(mass_matrix, sm_u.conj().T))
if no_bsm:
eg_vector = cardano_eqn(sm_ham)
else:
@@ -391,7 +438,7 @@ def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES,
scale_matrix = np.array(
[[0, 0, 0], [0, sc1, 0], [0, 0, sc2]]
)
- bsm_term = (energy**(dim-3)) * np.matmul(new_physics_u, np.matmul(scale_matrix, new_physics_u.conj().T))
+ bsm_term = (energy**(dim-3)) * np.dot(new_physics_u, np.dot(scale_matrix, new_physics_u.conj().T))
bsm_ham = sm_ham + bsm_term
eg_vector = cardano_eqn(bsm_ham)
diff --git a/utils/gf.py b/utils/gf.py
index c316ed2..79f1dd0 100644
--- a/utils/gf.py
+++ b/utils/gf.py
@@ -66,6 +66,7 @@ def steering_params(args):
# params.years = [999]
params.minFitEnergy = args.binning[0] # GeV
+ params.load_data_from_text_file = False
return params
diff --git a/utils/misc.py b/utils/misc.py
index 7e40d13..837c145 100644
--- a/utils/misc.py
+++ b/utils/misc.py
@@ -22,11 +22,6 @@ import numpy as np
from utils.enums import Likelihood
-DTYPE = np.float128
-CDTYPE = np.complex128
-PI = np.arccos(DTYPE(-1))
-
-
class SortingHelpFormatter(argparse.HelpFormatter):
"""Sort argparse help options alphabetically."""
def add_arguments(self, actions):
diff --git a/utils/multinest.py b/utils/multinest.py
index f3595f9..387ff90 100644
--- a/utils/multinest.py
+++ b/utils/multinest.py
@@ -95,6 +95,7 @@ def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter):
resume = False,
verbose = True
)
+ if result is None: return None
analyser = analyse.Analyzer(
outputfiles_basename=prefix, n_params=n_params
diff --git a/utils/plot.py b/utils/plot.py
index 65635bc..b63126b 100644
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -405,7 +405,7 @@ def plot_sens_fixed_angle(data, outfile, outformat, args):
print 'limit = {0}'.format(lim)
label = '{0} : {1} : {2}'.format(*misc_utils.solve_ratio(src))
if lim < yranges[0]: yranges[0] = lim
- if lim > yranges[1]: yranges[1] = lim+arr_len+1
+ if lim > yranges[1]: yranges[1] = lim+arr_len+2
# if lim > yranges[1]: yranges[1] = lim
xoff = 0.15
line = plt.Line2D(
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-17 15:10:23 +0000