show_pull_evolution.py
3.33 KB
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
#! /usr/bin/env python
# ==========================================================================
# This script displays the evolution of the pull distributions as function
# of the number of samples. The pull distributions are generated using the
# cspull script. The pull is defined by
#
# (fitted value - simulated value) / estimated error
#
# and the pull histogram allows to assess possible biases in the fit
# parameters and errors.
#
# Required 3rd party modules:
# - matplotlib
# - numpy
#
# Copyright (C) 2015 Juergen Knoedlseder
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# ==========================================================================
import matplotlib.pyplot as plt
import numpy as np
import sys
import csv
# ====================== #
# Read pull distribution #
# ====================== #
def read_pull(filename, parname):
"""
Read pull distribution.
"""
# Initialise list
values = []
# Open reader
reader = csv.reader(open(filename, 'r'), delimiter=',')
# Read rows
first = True
index = -1
for row in reader:
# Get column index if first row
if first:
try:
index = row.index(parname)
except:
sys.stdout.write('ERROR: Parameter "'+parname+'" not found in list:\n')
for p in row:
sys.stdout.write(' "'+p+'"\n')
raise NameError(parname)
# Handle data rows
else:
values.append(float(row[index]))
# Flag that first row has been passed
first = False
# Create numpy array
a = np.array(values)
# Return array
return a
# ======================== #
# Main routine entry point #
# ======================== #
if __name__ == '__main__':
"""
Display pull evolution generated by cspull using matplotlib.
"""
# Print usage information
usage = "Usage: show_pull_evolution filename parname"
if len(sys.argv) < 3:
sys.stdout.write(usage+"\n")
sys.exit()
# Extract parameters
filename = sys.argv[1]
parname = sys.argv[2]
# Read values from CSV file
values = read_pull(filename, parname)
# Compute mean and rms
nsamples = len(values)
samples = []
mean = []
rms = []
for i in range(1,nsamples):
samples.append(float(i+1))
mean.append(np.mean(values[0:i]))
rms.append(np.std(values[0:i]))
# Plot mean and rms evolution
plt.semilogx(samples, mean, 'r-', label="Mean")
plt.semilogx(samples, rms, 'b-', label="Std. deviation")
# Set plot
plt.xlabel("Samples")
plt.ylabel("Value")
plt.title(parname)
plt.grid(True)
# Set legend
plt.legend(loc="upper left")
# Show histogram
plt.show()