1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
|
#! /usr/bin/env python
# author : S. Mandalia
# s.p.mandalia@qmul.ac.uk
#
# date : April 10, 2018
"""
HESE BSM flavour ratio analysis script
"""
from __future__ import absolute_import, division
import numpy as np
import matplotlib as mpl
mpl.use('Agg')
from matplotlib import pyplot as plt
from matplotlib import rc
import GolemFitPy as gf
import fr
from utils import gf as gf_utils
from utils import likelihood as llh_utils
from utils import misc as misc_utils
from utils.enums import Likelihood, ParamTag
rc('text', usetex=False)
rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18})
RUN = True
z = 0+1E-6
scenarios = [
[np.sin(np.pi/2.)**2, z, z, z],
[z, np.cos(np.pi/2.)**4, z, z],
[z, z, np.sin(np.pi/2.)**2, z],
]
xticks = [r'$\mathcal{O}_{12}$', r'$\mathcal{O}_{13}$', r'$\mathcal{O}_{23}$']
def fit_flags(flag_dict):
flags = gf.FitParametersFlag()
for key in flag_dict.iterkeys():
flags.__setattr__(key, flag_dict[key])
return flags
default_flags = {
# False means it's not fixed in minimization
'astroFlavorAngle1' : True,
'astroFlavorAngle2' : True,
# 'astroENorm' : True,
# 'astroMuNorm' : True,
# 'astroTauNorm' : True,
'convNorm' : False,
'promptNorm' : False,
'muonNorm' : False,
'astroNorm' : False,
'astroParticleBalance' : True,
'astroDeltaGamma' : False,
'cutoffEnergy' : True,
'CRDeltaGamma' : False,
'piKRatio' : False,
'NeutrinoAntineutrinoRatio' : True,
'darkNorm' : True,
'domEfficiency' : True,
'holeiceForward' : True,
'anisotropyScale' : True,
'astroNormSec' : True,
'astroDeltaGammaSec' : True
}
def main():
args = fr.parse_args()
args.likelihood = Likelihood.GF_FREQ
fr.process_args(args)
misc_utils.print_args(args)
asimov_paramset, mcmc_paramset = fr.get_paramsets(args)
outfile = misc_utils.gen_outfile_name(args)
print '== {0:<25} = {1}'.format('outfile', outfile)
asimov_paramset = asimov_paramset.from_tag(ParamTag.BESTFIT)
mcmc_paramset = mcmc_paramset.from_tag(ParamTag.NUISANCE, invert=True)
sc_range = mcmc_paramset.from_tag(ParamTag.SCALE)[0].ranges
scan_scales = np.linspace(sc_range[0], sc_range[1], 100)
print 'scan_scales', scan_scales
if RUN:
datapaths = gf.DataPaths()
sparams = gf_utils.steering_params(args)
npp = gf.NewPhysicsParams()
fitter = gf.GolemFit(datapaths, sparams, npp)
fitter.SetFitParametersFlag(fit_flags(default_flags))
gf_utils.set_up_as(fitter, asimov_paramset)
x = []
y = []
mm_angles = mcmc_paramset.from_tag(ParamTag.MMANGLES)
for idx, scen in enumerate(scenarios):
scales = []
llhs = []
for yidx, an in enumerate(mm_angles):
an.value = scen[yidx]
for sc in scan_scales:
theta = scen + [sc]
print 'theta', theta
llh = llh_utils.triangle_llh(
theta=theta, args=args, asimov_paramset=asimov_paramset,
mcmc_paramset=mcmc_paramset, fitter=fitter
)
print 'llh', llh
scales.append(sc)
llhs.append(llh)
min_llh = np.min(llhs)
delta_llh = 2*(np.array(llhs) - min_llh)
print 'scales', scales
print 'delta_llh', delta_llh
n_arr = []
for i, d in enumerate(delta_llh):
# 90% for 1 DOF
if d < 2.71:
x.append(idx)
y.append(scales[i])
np.save('sens.npy', np.array([x, y]))
else:
x, y = np.load('sens.npy')
plt.plot(x, y, linewidth=0., marker='.')
plt.xticks(range(len(scenarios)), xticks)
plt.xlim(-1, len(scenarios))
plt.ylim(sc_range[0], sc_range[1])
plt.ylabel(r'${\rm log}_{10}\Lambda / GeV$')
plt.savefig('./sens.png', bbox_inches='tight', dpi=150)
main.__doc__ = __doc__
if __name__ == '__main__':
main()
|
