-
Notifications
You must be signed in to change notification settings - Fork 0
/
PyBLUE.py
executable file
·164 lines (124 loc) · 4.7 KB
/
PyBLUE.py
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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
#!/usr/bin/env python
import sys, os
from numpy import *
from ConfigParser import SafeConfigParser
############################################
def CalcCovariance( variation = 'm' ):
covtot = matrix( [ [0] * Nobservables ] * Nobservables )
for syst, rho in correlations.iteritems():
#if syst == "lumi": continue
#sv = array( [ unc[0][syst], unc[1][syst] ] )
slist = []
for channel in unc[variation]:
slist.append( channel[syst] )
sv = array( slist )
sm = diag( sv )
cov = sm * rho * sm
print "Source of uncertainty:",syst
print cov
covtot = covtot + cov
print "Total Covariance matrix for variation", variation, ":"
print covtot
return covtot
############################################
if len( sys.argv ) == 1:
print "Usage: ./PyBLUE config.blue"
exit(0)
parser = SafeConfigParser()
parser.read( sys.argv[1] )
Nobservables = int( parser.get( 'general', 'observables' ) )
print "INFO: Number of observables:", Nobservables
measurements_descriptions = parser.get( 'general', 'measurements' ).split()
Nmeasurements = len( measurements_descriptions )
print "INFO: Number of channels:", Nmeasurements
print "INFO: Channels:", measurements_descriptions
uncertainties_descriptions = parser.get( 'general', 'uncertainties' ).split()
Nuncertainties = len( uncertainties_descriptions )
print "INFO: Number of uncertainties:", Nuncertainties
print "INFO: Uncertainties:", uncertainties_descriptions
unc = {
'u' : [],
'm' : [],
'd' : []
}
blue_central = {}
blue_unc = {}
for variation in unc.keys():
blue_central[variation] = 0.
blue_unc[variation] = 0.
measurements = array( Nmeasurements * [ 0. ] )
n_c = 0
for channel in measurements_descriptions:
measurements[n_c] = float( parser.get( 'measurements', channel ) )
n_c += 1
unc_ch = {
'u' : {},
'm' : {},
'd' : {}
}
channel_u = channel + "_u"
channel_d = channel + "_d"
all_unc_this_ch_u = [ float(num) for num in parser.get( "uncertainties", channel_u ).split() ]
all_unc_this_ch_d = [ float(num) for num in parser.get( "uncertainties", channel_d ).split() ]
#all_unc_this_ch_m = [ 0.5 * ( up + down ) for up in all_unc_this_ch_u for down in all_unc_this_ch_d ]
all_unc_this_ch_m = []
for nx in range( len(all_unc_this_ch_u) ):
all_unc_this_ch_m.append( 0.5 * (all_unc_this_ch_u[nx] + all_unc_this_ch_d[nx]) )
n_u = 0
for s_unc in uncertainties_descriptions:
# print s_unc, n_u, channel, all_unc_this_ch_u[n_u]
unc_ch['u'][s_unc] = all_unc_this_ch_u[n_u]
unc_ch['m'][s_unc] = all_unc_this_ch_m[n_u]
unc_ch['d'][s_unc] = all_unc_this_ch_d[n_u]
n_u += 1
for variation in unc_ch.keys():
unc[variation].append( unc_ch[variation] )
print "INFO: Measurements:", measurements
print
correlations = {}
for s_unc in uncertainties_descriptions:
m = eval( parser.get( "correlations", s_unc ) )
correlations[s_unc] = matrix( m )
for variation in [ 'u', 'd', 'm' ]:
print "Variation", variation
cov = CalcCovariance( variation )
# pinv = pseudo-inverse. workaround for non-invertible matrices
#invcov = linalg.inv( cov )
invcov = linalg.pinv( cov )
print "Inverted covariance matrix:"
print invcov
unitvector = array( Nmeasurements * [ 1 ] )
l = dot( unitvector, invcov )
d = dot( unitvector, invcov )
d = dot( d, unitvector.transpose() )
# l = dot( unitvector.transpose(), invcov )
# d = dot( unitvector.transpose(), invcov )
# d = dot( d, unitvector )
l = l / d
print "Combination weights for variation", variation, ":"
print l
res = dot( l, measurements )
totunc = sqrt( dot( l, dot( cov, l.transpose() ) ) )
blue_central[variation] = float( res )
blue_unc[variation] = float( totunc )
print
print
print "/------------------------------------/"
print
print "Final results:"
print
for variation in [ 'u', 'd', 'm' ]:
print "Combination(%s) = %5.4f \pm %10.7f" % ( variation, blue_central[variation], blue_unc[variation] )
print
# AIB
#R_u = blue_unc['u'] / ( blue_unc['u'] + blue_unc['d'] )
#sigma_u = 2. * R_u * blue_unc['m']
#sigma_d = 2. * ( 1. - R_u ) * blue_unc['m']
#
#print "Combination(AIB) = %5.4f +%5.4f -%5.4f" % ( blue_central['m'], sigma_u, sigma_d )
denom = 0.5 * ( blue_unc['u'] + blue_unc['d'] )
factor = blue_unc['m'] / denom
upperError = factor * blue_unc['u']
lowerError = factor * blue_unc['d']
print "Combination(AIB) = %10.7f +%10.7f -%10.7f" % ( blue_central['m'], upperError, lowerError )
#print "Combination(AIB) = %10.7f +%8.7f -%8.7f" % ( blue_central['m'], blue_unc['u'], blue_unc['d'] )