aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rwxr-xr-xplot_sens.py29
-rwxr-xr-xsens.py5
-rw-r--r--submitter/mcmc_dag.py5
-rw-r--r--submitter/mcmc_submit.sub3
-rw-r--r--submitter/sens_dag.py29
-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
10 files changed, 62 insertions, 73 deletions
diff --git a/plot_sens.py b/plot_sens.py
index 2e3b1c6..9aa5089 100755
--- a/plot_sens.py
+++ b/plot_sens.py
@@ -145,10 +145,6 @@ def parse_args(args=None):
help='Make sensitivity plots'
)
parser.add_argument(
- '--data', type=misc_utils.parse_bool, default='False',
- help='Use HESE7 data'
- )
- parser.add_argument(
'--plot-statistic', type=misc_utils.parse_bool, default='False',
help='Plot MultiNest evidence or LLH value'
)
@@ -231,13 +227,11 @@ def main():
infile += '/golemfit/'
elif args.likelihood is Likelihood.GAUSSIAN:
infile += '/gaussian/'
- # infile += '/DIM{0}/fix_ifr/HESESim'.format(dim)
- # infile += '/DIM{0}/fix_ifr/'.format(dim)
- infile += '/DIM{0}/fix_ifr/data'.format(dim)
if args.likelihood is Likelihood.GAUSSIAN:
infile += '{0}/'.format(str(args.sigma_ratio).replace('.', '_'))
- infile += 'fr_stat/{0}/{1}/fr_stat'.format(args.stat_method, args.run_method) \
- + misc_utils.gen_identifier(argsc)
+ infile += '/{0}/{1}/{2}/fr_stat'.format(
+ *map(misc_utils.parse_enum, [args.stat_method, args.run_method, args.data])
+ ) + misc_utils.gen_identifier(argsc)
print '== {0:<25} = {1}'.format('infile', infile)
if args.split_jobs:
@@ -279,17 +273,14 @@ def main():
argsc.dimension = dim
_, scale_region = fr_utils.estimate_scale(argsc)
argsc.scale_region = scale_region
- # infile = base_infile + '/DIM{0}/fix_ifrHESESim/'.format(dim)
- # infile = base_infile + '/DIM{0}/fix_ifr/'.format(dim)
- infile = base_infile + '/DIM{0}/fix_ifr/data'.format(dim)
if args.likelihood is Likelihood.GAUSSIAN:
infile += '{0}/'.format(str(args.sigma_ratio).replace('.', '_'))
- infile += 'fr_stat/'
for isrc, src in enumerate(args.source_ratios):
argsc.source_ratio = src
- finfile = infile +'/{0}/{1}/fr_stat'.format(args.stat_method, args.run_method) \
- + misc_utils.gen_identifier(argsc)
+ finfile = infile +'/{0}/{1}/{2}/fr_stat'.format(
+ *map(misc_utils.parse_enum, [args.stat_method, args.run_method, args.data])
+ ) + misc_utils.gen_identifier(argsc)
basename = os.path.dirname(finfile) + '/statrun/'
baseoutfile = basename[:5]+basename[5:].replace('data', 'plots')
if args.data:
@@ -343,11 +334,9 @@ def main():
print 'Plotting sensitivities'
basename = args.infile[:5]+args.infile[5:].replace('data', 'plots')
- if args.likelihood in [Likelihood.GOLEMFIT, Likelihood.GF_FREQ]:
- baseoutfile = basename + 'golemfit/{1}'.format(args.likelihood, args.stat_method)
- if args.data: baseoutfile += '/data'
- elif args.likelihood is Likelihood.GAUSSIAN:
- baseoutfile = basename + 'gaussian/{1}'.format(args.likelihood, args.stat_method)
+ baseoutfile = basename + '/{1}/{2}/{3}'.format(
+ *map(misc_utils.parse_enum, [args.likelihood, args.stat_method, args.data])
+ )
if args.run_method is SensitivityCateg.FULL:
plot_utils.plot_sens_full(
diff --git a/sens.py b/sens.py
index 6145f4a..581a3ca 100755
--- a/sens.py
+++ b/sens.py
@@ -207,8 +207,9 @@ def main():
eval_dim = args.sens_bins
else: eval_dim = 1
- out = args.outfile+'/{0}/{1}/fr_stat'.format(args.stat_method, args.run_method) \
- + misc_utils.gen_identifier(args)
+ out = args.outfile+'/{0}/{1}/{2}/fr_stat'.format(
+ *map(misc_utils.parse_enum, [args.stat_method, args.run_method, args.data])
+ ) + misc_utils.gen_identifier(args)
if args.sens_run:
if args.likelihood in [Likelihood.GOLEMFIT, Likelihood.GF_FREQ]:
fitter = gf_utils.setup_fitter(args, asimov_paramset)
diff --git a/submitter/mcmc_dag.py b/submitter/mcmc_dag.py
index 7a0cde9..e40c043 100644
--- a/submitter/mcmc_dag.py
+++ b/submitter/mcmc_dag.py
@@ -38,7 +38,7 @@ dimension = [3]
# dimension = [4, 5, 7, 8]
GLOBAL_PARAMS.update(dict(
threads = 1,
- binning = '6e4 1e7 5',
+ binning = '6e4 1e7 20',
no_bsm = 'False',
scale_region = "1E10",
energy_dependance = 'spectral',
@@ -57,7 +57,7 @@ GLOBAL_PARAMS.update(dict(
# GolemFit
GLOBAL_PARAMS.update(dict(
ast = 'p2_0',
- data = 'real'
+ data = 'asimov'
))
# Plot
@@ -82,6 +82,7 @@ with open(outfile, 'w') as f:
outchains = outchain_head + '/fix_ifr/'
if GLOBAL_PARAMS['likelihood'].lower() == 'gaussian':
outchains += '{0}/'.format(str(GLOBAL_PARAMS['sigma_ratio']).replace('.', '_'))
+ outchains += '{0}/'.format(GLOBAL_PARAMS['data'].lower())
outchains += 'mcmc_chain'
f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
f.write('VARS\tjob{0}\tdimension="{1}"\n'.format(job_number, dim))
diff --git a/submitter/mcmc_submit.sub b/submitter/mcmc_submit.sub
index 998932e..fa36250 100644
--- a/submitter/mcmc_submit.sub
+++ b/submitter/mcmc_submit.sub
@@ -32,9 +32,10 @@ Notification = never
# +TransferOutput=""
++NATIVE_OS = True
# Requirements = IS_GLIDEIN && HAS_CVMFS_icecube_opensciencegrid_org && (OpSysAndVer =?= "CentOS6" || OpSysAndVer =?= "RedHat6" || OpSysAndVer =?= "SL6")
# Requirements = IS_GLIDEIN
-Requirements = (OpSysMajorVer =?= 6)
+# Requirements = (OpSysMajorVer =?= 6)
# GO!
queue
diff --git a/submitter/sens_dag.py b/submitter/sens_dag.py
index 666a0dd..5e00dfd 100644
--- a/submitter/sens_dag.py
+++ b/submitter/sens_dag.py
@@ -35,27 +35,25 @@ GLOBAL_PARAMS.update(dict(
# MultiNest
GLOBAL_PARAMS.update(dict(
- mn_live_points = 300,
+ mn_live_points = 1000,
mn_tolerance = 0.1,
mn_output = './mnrun'
))
# FR
-# dimension = [3]
-dimension = [3, 6]
-# dimension = [4, 5, 7, 8]
-# dimension = [3, 4, 5, 6, 7, 8]
+dimension = [3]
+# dimension = [3, 6]
GLOBAL_PARAMS.update(dict(
threads = 1,
- # binning = '6e4 1e7 5',
- binning = '1e5 1e7 5',
+ binning = '6e4 1e7 20',
+ # binning = '1e5 1e7 5',
no_bsm = 'False',
scale_region = "1E10",
energy_dependance = 'spectral',
spectral_index = -2,
fix_mixing = 'False',
fix_mixing_almost = 'False',
- fold_index = 'False'
+ fold_index = 'True'
))
# Likelihood
@@ -67,7 +65,7 @@ GLOBAL_PARAMS.update(dict(
# GolemFit
GLOBAL_PARAMS.update(dict(
ast = 'p2_0',
- data = 'real'
+ data = 'asimov'
))
# Plot
@@ -75,8 +73,9 @@ GLOBAL_PARAMS.update(dict(
plot_statistic = 'True'
))
-outfile = 'dagman_FR_SENS_{0}_{1}_{2}_data_100TeV.submit'.format(
- GLOBAL_PARAMS['stat_method'], GLOBAL_PARAMS['run_method'], GLOBAL_PARAMS['likelihood']
+outfile = 'dagman_FR_SENS_{0}_{1}_{2}_{3}.submit'.format(
+ GLOBAL_PARAMS['stat_method'], GLOBAL_PARAMS['run_method'],
+ GLOBAL_PARAMS['likelihood'], GLOBAL_PARAMS['data']
)
golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit'
condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/sens_submit.sub'
@@ -97,14 +96,9 @@ with open(outfile, 'w') as f:
)
for frs in fix_sfr_mfr:
print 'frs', frs
- # output = outchain_head + '/fix_ifr/'
- # output = outchain_head + '/fix_ifr/HESESim'
- # output = outchain_head + '/fix_ifr/sim'
- # output = outchain_head + '/fix_ifr/data'
- output = outchain_head + '/fix_ifr/data/100TeV/'
+ output = outchain_head + '/fix_ifr/'
if GLOBAL_PARAMS['likelihood'].lower() == 'gaussian':
output += '{0}/'.format(str(GLOBAL_PARAMS['sigma_ratio']).replace('.', '_'))
- output += 'fr_stat'
for r in xrange(sens_runs):
print 'run', r
f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
@@ -130,7 +124,6 @@ with open(outfile, 'w') as f:
output = outchain_head + '/full/'
if GLOBAL_PARAMS['likelihood'].lower() == 'gaussian':
output += '{0}/'.format(str(GLOBAL_PARAMS['sigma_ratio']).replace('.', '_'))
- output += 'fr_stat'
for r in xrange(sens_runs):
print 'run', r
f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
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 09:49:31 +0000