forked from maikeJung/LLH_M2
-
Notifications
You must be signed in to change notification settings - Fork 0
/
PLOT_HIST.py
119 lines (98 loc) · 4.62 KB
/
PLOT_HIST.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
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d
import matplotlib.patches as mpatches
### PSEUDO EXPERIMENTS ###
# read in data
def readData(filename, input_mass, bin_width):
h = open(filename,'r')
masses = []
for line in h:
m = float(line)
masses.append(m)
h.close()
events = len(masses)
# store values in histogram -> values_pseudo: values of bins; mass_hist_pseudo: center of the bins
bins = np.arange(input_mass-2.0 - bin_width/2.0, input_mass+2.0, bin_width)
values_pseudo, m = np.histogram(masses, bins=bins)
mass_hist_pseudo = np.arange(input_mass-2.0, input_mass+2.0 - bin_width/2.0, bin_width)
return mass_hist_pseudo, values_pseudo/float(events)
#return mass_hist_pseudo, values_pseudo
def calcError(mass, events, distance, bin_width):
# calculate the 1 sigma uncertainty, by checking where 68% of events are detected
mass_hist, values = readData("DATA/masses2_"+str(distance)+"Mpc_"+str(events)+".0Events_"+str(mass)+"eV2.txt", mass, bin_width)
mass_hist = mass_hist[::-1]
values = values[::-1]
values = np.cumsum(values)
f = interp1d(values, mass_hist)
#return f(0.84135), f(0.5), f(0.15865)
# 90%C.L., 2sigma, 1sigma and mean
return f(0.97725),f(0.9), f(0.84135), f(0.5), f(0.15865), f(0.02275)
def plotHist(mass, events, distance, bin_width):
mass_hist, values = readData("DATA/masses2_"+str(distance)+"Mpc_"+str(events)+".0Events_"+str(mass)+"eV2.txt", mass, bin_width)
# create a histogram of pseudo experiments
fig=plt.figure()
ax = fig.add_subplot(111)
plt.gcf().subplots_adjust(top=0.85) # to be able to read the lables
plt.gcf().subplots_adjust(bottom=0.15)
plt.gcf().subplots_adjust(left=0.15)
axes.yaxis.set_tick_params(labelsize=17, width=1)
plt.bar(mass_hist-bin_width/2.0, values, bin_width, color="blue", label=str(distance)+"Mpc "+str(events)+"Events \n"+str(mass)+"eV", alpha=0.6)
#plt.xlim(0.6,1.4)
plt.xlabel('mass [eV]', fontsize=18)
plt.ylabel('# of pseudo experiments', fontsize=18)
plt.legend(loc='upper right', fontsize=14)
#fig.savefig('PLOTS10000/'+str(mass)+'eV_'+str(distance)+'Mpc_'+str(events)+'Events.png')
def calculate_error(input_mass, events, distance):
lower = []
value = []
upper = []
lower2sigma = []
upper2sigma = []
low90CL = []
for mass in input_mass:
low2sigma, t90CL, low, val, upp, upp2sigma = calcError(mass, int(events), distance, 0.0001)
lower.append(low)
value.append(val)
upper.append(upp)
lower2sigma.append(low2sigma)
upper2sigma.append(upp2sigma)
low90CL.append(t90CL)
lower = np.asarray(lower)
value = np.asarray(value)
upper = np.asarray(upper)
lower2sigma = np.asarray(lower2sigma)
upper2sigma = np.asarray(upper2sigma)
low90CL = np.asarray(low90CL)
return lower2sigma, low90CL, lower, value, upper, upper2sigma
#set parameters
distance = 1.0
events = 160
input_mass = [0.0, 0.01, 0.3, 0.5]
lower2sigma, test90CL, lower, value, upper, upper2sigma = calculate_error(input_mass, events, distance)
#plot
fig=plt.figure()
plt.gcf().subplots_adjust(top=0.85) # to be able to read the lables
plt.gcf().subplots_adjust(bottom=0.15)
plt.gcf().subplots_adjust(left=0.15)
minput, = plt.plot(input_mass,input_mass,color='#000000', linewidth=1.5, linestyle='--', label=r'$m^2_{\mathrm{input}}=m^2_{\mathrm{found}}$')
mfound, = plt.plot(input_mass,value,color='red', linewidth=1.5, linestyle='--', label=r'$<m^2_{\mathrm{found}}>$')
m90CL, = plt.plot(input_mass,test90CL,color='blue', linewidth=1.5, linestyle='-', label=r'$90\%$ C.L.')
green_patch = mpatches.Patch(color='#66ff33', label=r'$1\,\sigma$')
yellow_patch = mpatches.Patch(color='#ffff00', label=r'$2\,\sigma$')
plt.fill_between(input_mass, lower2sigma, upper2sigma,edgecolor='#ffff00', facecolor='#ffff00')
plt.fill_between(input_mass, lower, upper,edgecolor='#66ff33', facecolor='#66ff33')
plt.title('Uncertainty of determined neutrino mass \n' + r'real detector with constant noise floor', fontsize=16)
plt.xlabel(r'$m^2_{\mathrm{input}}$ [eV$^2$]', fontsize=16)
plt.ylabel(r'reconstructed mass $m^2_{\mathrm{found}}$ [eV$^2$]', fontsize=16)
plt.legend(loc='upper left', numpoints = 1, fontsize=14.5, handles=[green_patch, yellow_patch, minput, mfound, m90CL])
axes = plt.gca()
axes.set_xlim([0.0,1.501])
axes.set_ylim([0.0,1.6])
axes.xaxis.set_tick_params(labelsize=15, width=1)
axes.yaxis.set_tick_params(labelsize=15, width=1)
#plot histogram for certain configuration
hit_width = 0.06
plotHist(0.3,events,distance, hit_width)
plotHist(0.5,events,distance, hit_width)
plt.show()