From 729ebbcc3da2a194fee324b5b0108a27a2a8c091 Mon Sep 17 00:00:00 2001
From: knodlseder <>
Date: Wed, 3 Aug 2011 23:33:06 +0000
Subject: [PATCH] initial version

---
 examples/show_ts_distribution.py | 155 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 155 insertions(+), 0 deletions(-)
 create mode 100755 examples/show_ts_distribution.py

diff --git a/examples/show_ts_distribution.py b/examples/show_ts_distribution.py
new file mode 100755
index 0000000..4f4a40d
--- /dev/null
+++ b/examples/show_ts_distribution.py
@@ -0,0 +1,155 @@
+#! /usr/bin/env python
+# ==========================================================================
+# This script displays pull histograms that have been 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) 2011 Jurgen Knodlseder
+#
+# 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.mlab as mlab
+import matplotlib.pyplot as plt
+import numpy as np
+import sys
+import csv
+import math
+
+
+# ==================== #
+# Read TS distribution #
+# ==================== #
+def read_ts(filename):
+	"""
+	Read TS 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("TS")
+			except:
+				print 'ERROR: Column "TS" not found in file'
+				raise NameError("TS")
+			#print index
+
+		# 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)
+	#print a
+	
+	# Return array
+	return a
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+	"""
+	Display pull histogram generated by cspull using matplotlib.
+	"""
+	# Print usage information
+	usage = "Usage: show_ts_distribution filename [bins]"
+	if len(sys.argv) < 2:
+		print usage
+		sys.exit()
+
+	# Extract parameters
+	filename = sys.argv[1]
+	if len(sys.argv) == 3:
+		nbins = int(sys.argv[2])
+	else:
+		nbins = 30
+	
+	# Set lin/log
+	log = False
+	
+	# Read values from CSV file
+	values = read_ts(filename)
+	print len(values), "values read."
+	
+	# Set range and adapt number of bins. We make sure that we have a
+	# histogram bin centred on 0 that will capture the TS=0 part.
+	max_value = max(values)
+	binsize   = max_value/nbins
+	min_value = 0.5*binsize
+	while min_value > min(values):
+		min_value -= binsize
+		nbins     += 1
+	
+	# Create histogram
+	n, bins, patches = plt.hist(values, nbins, range=[min_value, max_value], \
+	                            align='mid', normed=True, facecolor='green', \
+								log=log)
+
+	# Create expected distribution (onyl for positive TS). We compute the
+	# full and half of the distribution, as for positively constrained
+	# quantities half of the values have TS=0
+	x      = []
+	y_full = []
+	y_half = []
+	norm   = 1.0/(math.sqrt(2.0)*math.sqrt(math.pi))
+	width  = bins[1]-bins[0]
+	x.append(0.0)
+	y_full.append(1.0/width)
+	y_half.append(0.5/width)
+	for bin in bins:
+		ts = bin + 0.5*width
+		if ts > 0:
+			y  = norm*math.pow(ts,-0.5)*math.exp(-0.5*ts)
+			x.append(ts)
+			y_full.append(y)
+			y_half.append(0.5*y)
+	
+	# Show expected distribution (half)
+	if log:
+		l = plt.semilogy(x, y_half, 'ro', linewidth=2)
+	else:
+		l = plt.plot(x, y_half, 'r-', linewidth=2)
+
+	# Set plot
+	plt.xlabel('TS')
+	plt.ylabel('Arbitrary units')
+	plt.title('TS distribution')
+	plt.grid(True)
+
+	# Show histogram
+	plt.show()
+	
\ No newline at end of file
--
libgit2 0.21.4