Andreas Specovius / ctools

Commit 729ebbcc3da2a194fee324b5b0108a27a2a8c091

Authored by knodlseder
1 parent 130f3114

initial version

Inline Side-by-side
Showing 1 changed file with 155 additions and 0 deletions
examples/show_ts_distribution.py 0 → 100755
  View file @729ebbc
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script displays pull histograms that have been generated using the
  4 +# cspull script. The pull is defined by
  5 +#
  6 +# (fitted value - simulated value) / estimated error
  7 +#
  8 +# and the pull histogram allows to assess possible biases in the fit
  9 +# parameters and errors.
  10 +#
  11 +# Required 3rd party modules:
  12 +# - matplotlib
  13 +# - numpy
  14 +#
  15 +# Copyright (C) 2011 Jurgen Knodlseder
  16 +#
  17 +# This program is free software: you can redistribute it and/or modify
  18 +# it under the terms of the GNU General Public License as published by
  19 +# the Free Software Foundation, either version 3 of the License, or
  20 +# (at your option) any later version.
  21 +#
  22 +# This program is distributed in the hope that it will be useful,
  23 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  24 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  25 +# GNU General Public License for more details.
  26 +#
  27 +# You should have received a copy of the GNU General Public License
  28 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  29 +#
  30 +# ==========================================================================
  31 +import matplotlib.mlab as mlab
  32 +import matplotlib.pyplot as plt
  33 +import numpy as np
  34 +import sys
  35 +import csv
  36 +import math
  37 +
  38 +
  39 +# ==================== #
  40 +# Read TS distribution #
  41 +# ==================== #
  42 +def read_ts(filename):
  43 + """
  44 + Read TS distribution.
  45 + """
  46 + # Initialise list
  47 + values = []
  48 +
  49 + # Open reader
  50 + reader = csv.reader(open(filename, 'r'), delimiter=',')
  51 +
  52 + # Read rows
  53 + first = True
  54 + index = -1
  55 + for row in reader:
  56 +
  57 + # Get column index if first row
  58 + if first:
  59 + try:
  60 + index = row.index("TS")
  61 + except:
  62 + print 'ERROR: Column "TS" not found in file'
  63 + raise NameError("TS")
  64 + #print index
  65 +
  66 + # Handle data rows
  67 + else:
  68 + values.append(float(row[index]))
  69 +
  70 + # Flag that first row has been passed
  71 + first = False
  72 +
  73 + # Create numpy array
  74 + a = np.array(values)
  75 + #print a
  76 +
  77 + # Return array
  78 + return a
  79 +
  80 +
  81 +# ======================== #
  82 +# Main routine entry point #
  83 +# ======================== #
  84 +if __name__ == '__main__':
  85 + """
  86 + Display pull histogram generated by cspull using matplotlib.
  87 + """
  88 + # Print usage information
  89 + usage = "Usage: show_ts_distribution filename [bins]"
  90 + if len(sys.argv) < 2:
  91 + print usage
  92 + sys.exit()
  93 +
  94 + # Extract parameters
  95 + filename = sys.argv[1]
  96 + if len(sys.argv) == 3:
  97 + nbins = int(sys.argv[2])
  98 + else:
  99 + nbins = 30
  100 +
  101 + # Set lin/log
  102 + log = False
  103 +
  104 + # Read values from CSV file
  105 + values = read_ts(filename)
  106 + print len(values), "values read."
  107 +
  108 + # Set range and adapt number of bins. We make sure that we have a
  109 + # histogram bin centred on 0 that will capture the TS=0 part.
  110 + max_value = max(values)
  111 + binsize = max_value/nbins
  112 + min_value = 0.5*binsize
  113 + while min_value > min(values):
  114 + min_value -= binsize
  115 + nbins += 1
  116 +
  117 + # Create histogram
  118 + n, bins, patches = plt.hist(values, nbins, range=[min_value, max_value], \
  119 + align='mid', normed=True, facecolor='green', \
  120 + log=log)
  121 +
  122 + # Create expected distribution (onyl for positive TS). We compute the
  123 + # full and half of the distribution, as for positively constrained
  124 + # quantities half of the values have TS=0
  125 + x = []
  126 + y_full = []
  127 + y_half = []
  128 + norm = 1.0/(math.sqrt(2.0)*math.sqrt(math.pi))
  129 + width = bins[1]-bins[0]
  130 + x.append(0.0)
  131 + y_full.append(1.0/width)
  132 + y_half.append(0.5/width)
  133 + for bin in bins:
  134 + ts = bin + 0.5*width
  135 + if ts > 0:
  136 + y = norm*math.pow(ts,-0.5)*math.exp(-0.5*ts)
  137 + x.append(ts)
  138 + y_full.append(y)
  139 + y_half.append(0.5*y)
  140 +
  141 + # Show expected distribution (half)
  142 + if log:
  143 + l = plt.semilogy(x, y_half, 'ro', linewidth=2)
  144 + else:
  145 + l = plt.plot(x, y_half, 'r-', linewidth=2)
  146 +
  147 + # Set plot
  148 + plt.xlabel('TS')
  149 + plt.ylabel('Arbitrary units')
  150 + plt.title('TS distribution')
  151 + plt.grid(True)
  152 +
  153 + # Show histogram
  154 + plt.show()
  155 +
0 156 \ No newline at end of file
... ...