diff options
Diffstat (limited to 'mcmc_scan.py')
| -rwxr-xr-x | mcmc_scan.py | 24 |
1 files changed, 22 insertions, 2 deletions
diff --git a/mcmc_scan.py b/mcmc_scan.py index 2183280..688345e 100755 --- a/mcmc_scan.py +++ b/mcmc_scan.py @@ -29,6 +29,7 @@ FIX_MIXING = False FIX_SFR = True SOURCE_FR = [1, 2, 0] FIX_SCALE = True +NO_BSM = False DIMENSION = 3 ENERGY = 1000000 # GeV @@ -80,7 +81,7 @@ def set_up_as(golemfit, parms): golemfit.SetupAsimov(asimov_parameters) -def get_best_fit(golemfit, parms): +def get_llh(golemfit, parms): fitparams = gf.FitParameters() fitparams.convNorm = parms['convNorm'] fitparams.promptNorm = parms['promptNorm'] @@ -202,6 +203,12 @@ def params_to_BSMu(theta): [[0, 0, 0], [0, MASS_EIGENVALUES[0], 0], [0, 0, MASS_EIGENVALUES[1]]] ) sm_ham = (1./(2*ENERGY))*np.dot(NUFIT_U, np.dot(mass_matrix, NUFIT_U.conj().T)) + if NO_BSM: + eg_vector = cardano_eqn(sm_ham) + tu = test_uni(eg_vector) + if not abs(np.trace(tu) - 3.) < 1e-5 or not abs(np.sum(tu) - 3.) < 1e-5: + raise AssertionError('Matrix is not unitary!\neg_vector\n{0}\ntest u\n{1}'.format(eg_vector, tu)) + return eg_vector new_physics_u = angles_to_u((s12_2, c13_4, s23_2, dcp)) scale_matrix = np.array( @@ -255,7 +262,7 @@ def triangle_llh(theta): if FLAT: return 10. else: - return get_best_fit(GOLEMFIT, fit_values) + return get_llh(GOLEMFIT, fit_values) def lnprior(theta): @@ -336,6 +343,10 @@ def parse_args(): help='Set the source flavour ratio for the case when you want to fix it' ) parser.add_argument( + '--no-bsm', type=str, default='False', + help='Turn off BSM terms' + ) + parser.add_argument( '--fix-mixing', type=str, default='False', help='Fix all mixing parameters except one' ) @@ -399,6 +410,7 @@ def main(): global FIX_SCALE global SCALE global SCALE2_BOUNDS + global NO_BSM DIMENSION = args.dimension ENERGY = args.energy @@ -434,6 +446,13 @@ def main(): else: raise ValueError + if args.no_bsm.lower() == 'true': + NO_BSM = True + elif args.no_bsm.lower() == 'false': + NO_BSM = False + else: + raise ValueError + if FIX_SFR: SOURCE_FR = np.array(args.source_ratio) / float(np.sum(args.source_ratio)) @@ -456,6 +475,7 @@ def main(): print 'MEASURED_FR = {0}'.format(MEASURED_FR) print 'SIGMA = {0}'.format(SIGMA) print 'FLAT = {0}'.format(FLAT) + print 'NO_BSM = {0}'.format(NO_BSM) print 'ENERGY = {0}'.format(ENERGY) print 'DIMENSION = {0}'.format(DIMENSION) print 'SCALE2_BOUNDS = {0}'.format(SCALE2_BOUNDS) |