aboutsummaryrefslogtreecommitdiffstats
path: root/utils
diff options
context:
space:
mode:
Diffstat (limited to 'utils')
-rw-r--r--utils/enums.py5
-rw-r--r--utils/fr.py39
-rw-r--r--utils/gf.py14
-rw-r--r--utils/likelihood.py2
-rw-r--r--utils/misc.py4
5 files changed, 34 insertions, 30 deletions
diff --git a/utils/enums.py b/utils/enums.py
index b80b712..c40d681 100644
--- a/utils/enums.py
+++ b/utils/enums.py
@@ -11,9 +11,8 @@ from enum import Enum
class DataType(Enum):
- REAL = 1
- FAKE = 2
- ASMIMOV = 3
+ REAL = 1
+ ASIMOV = 2
class EnergyDependance(Enum):
diff --git a/utils/fr.py b/utils/fr.py
index 64e2706..639927f 100644
--- a/utils/fr.py
+++ b/utils/fr.py
@@ -140,16 +140,16 @@ def cardano_eqn(ham):
)
a = -np.trace(ham)
- b = (0.5) * ((np.trace(ham))**2 - np.trace(np.dot(ham, ham)))
+ b = DTYPE(1)/2 * ((np.trace(ham))**DTYPE(2) - np.trace(np.dot(ham, ham)))
c = -np.linalg.det(ham)
- Q = (1/9.) * (a**2 - 3*b)
- R = (1/54.) * (2*a**3 - 9*a*b + 27*c)
- theta = np.arccos(R / np.sqrt(Q**3))
+ 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)))
- E1 = -2 * np.sqrt(Q) * np.cos(theta/3.) - (1/3.)*a
- E2 = -2 * np.sqrt(Q) * np.cos((theta - 2.*PI)/3.) - (1/3.)*a
- E3 = -2 * np.sqrt(Q) * np.cos((theta + 2.*PI)/3.) - (1/3.)*a
+ 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
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]
@@ -383,7 +383,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.dot(sm_u, np.dot(mass_matrix, sm_u.conj().T))
+ sm_ham = (1./(2*energy))*np.matmul(sm_u, np.matmul(mass_matrix, sm_u.conj().T))
if no_bsm:
eg_vector = cardano_eqn(sm_ham)
else:
@@ -391,23 +391,16 @@ 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.dot(new_physics_u, np.dot(scale_matrix, new_physics_u.conj().T))
-
+ bsm_term = (energy**(dim-3)) * np.matmul(new_physics_u, np.matmul(scale_matrix, new_physics_u.conj().T))
bsm_ham = sm_ham + bsm_term
eg_vector = cardano_eqn(bsm_ham)
if check_uni:
- tu = test_unitarity(eg_vector)
- if not np.abs(np.trace(tu) - 3.) < epsilon or \
- not np.abs(np.sum(tu) - 3.) < epsilon:
- raise AssertionError(
- 'Matrix is not unitary!\neg_vector\n{0}\ntest '
- 'u\n{1}'.format(eg_vector, tu)
- )
+ test_unitarity(eg_vector, rse=True, epsilon=epsilon)
return eg_vector
-def test_unitarity(x, prnt=False):
+def test_unitarity(x, prnt=False, rse=False, epsilon=None):
"""Test the unitarity of a matrix.
Parameters
@@ -418,6 +411,9 @@ def test_unitarity(x, prnt=False):
prnt : bool
Print the result
+ rse : bool
+ Raise Assertion if matrix is not unitary
+
Returns
----------
numpy ndarray
@@ -435,6 +431,13 @@ def test_unitarity(x, prnt=False):
f = np.abs(np.dot(x, x.conj().T), dtype=DTYPE)
if prnt:
print 'Unitarity test:\n{0}'.format(f)
+ if rse:
+ if not np.abs(np.trace(f) - 3.) < epsilon or \
+ not np.abs(np.sum(f) - 3.) < epsilon:
+ raise AssertionError(
+ 'Matrix is not unitary!\nx\n{0}\ntest '
+ 'u\n{1}'.format(x, f)
+ )
return f
diff --git a/utils/gf.py b/utils/gf.py
index 7ded152..cc093e1 100644
--- a/utils/gf.py
+++ b/utils/gf.py
@@ -55,7 +55,6 @@ def fit_flags(llh_paramset):
def steering_params(args):
steering_categ = args.ast
- # params = gf.SteeringParams(gf.sampleTag.HESE)
params = gf.SteeringParams(gf.sampleTag.MagicTau)
params.quiet = False
params.fastmode = True
@@ -66,14 +65,13 @@ def steering_params(args):
# For Tianlu
# params.years = [999]
- params.minFitEnergy = 1.0e5 # GeV
+ # params.minFitEnergy = 1.0e5 # GeV
return params
def set_up_as(fitter, params):
print 'Injecting the model', params
- # asimov_params = gf.FitParameters(gf.sampleTag.HESE)
asimov_params = gf.FitParameters(gf.sampleTag.MagicTau)
for parm in params:
asimov_params.__setattr__(parm.name, float(parm.value))
@@ -85,13 +83,14 @@ def setup_fitter(args, asimov_paramset):
sparams = steering_params(args)
npp = gf.NewPhysicsParams()
fitter = gf.GolemFit(datapaths, sparams, npp)
- # comment to use data
- # set_up_as(fitter, asimov_paramset)
+ if args.data is DataType.ASIMOV:
+ set_up_as(fitter, asimov_paramset)
+ elif args.data is DataType.REAL:
+ print 'Using MagicTau DATA'
return fitter
def get_llh(fitter, params):
- # fitparams = gf.FitParameters(gf.sampleTag.HESE)
fitparams = gf.FitParameters(gf.sampleTag.MagicTau)
for parm in params:
fitparams.__setattr__(parm.name, float(parm.value))
@@ -101,7 +100,6 @@ def get_llh(fitter, params):
def get_llh_freq(fitter, params):
print 'setting to {0}'.format(params)
- # fitparams = gf.FitParameters(gf.sampleTag.HESE)
fitparams = gf.FitParameters(gf.sampleTag.MagicTau)
for parm in params:
fitparams.__setattr__(parm.name, float(parm.value))
@@ -118,7 +116,7 @@ def data_distributions(fitter):
def gf_argparse(parser):
parser.add_argument(
- '--data', default='real', type=partial(enum_parse, c=DataType),
+ '--data', default='asimov', type=partial(enum_parse, c=DataType),
choices=DataType, help='select datatype'
)
parser.add_argument(
diff --git a/utils/likelihood.py b/utils/likelihood.py
index ebc3539..ed598c8 100644
--- a/utils/likelihood.py
+++ b/utils/likelihood.py
@@ -146,7 +146,7 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter):
u = fr_utils.params_to_BSMu(
theta = bsm_angles,
dim = args.dimension,
- energy = args.energy,
+ energy = bin_centers[i_sf],
mass_eigenvalues = mass_eigenvalues,
sm_u = sm_u,
no_bsm = args.no_bsm,
diff --git a/utils/misc.py b/utils/misc.py
index 5483675..7e40d13 100644
--- a/utils/misc.py
+++ b/utils/misc.py
@@ -106,6 +106,10 @@ def parse_bool(s):
raise ValueError
+def parse_enum(e):
+ return '{0}'.format(e).split('.')[1].lower()
+
+
def print_args(args):
"""Print the input arguments.
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-17 18:40:30 +0000