Lili Yang / ctools

Commit 055de2c798d2ff229a7ddfed5de03a2d32b8fea5

Authored by knodlseder
1 parent dd4efb46

move here from test directory

Inline Side-by-side
Showing 1 changed file with 125 additions and 0 deletions
examples/pull_histogram.py 0 → 100755
  View file @055de2c
  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 +# Requirements:
  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 +
  37 +
  38 +# ====================== #
  39 +# Read pull distribution #
  40 +# ====================== #
  41 +def read_pull(filename, parname):
  42 + """
  43 + Read pull distribution.
  44 + """
  45 + # Initialise list
  46 + values = []
  47 +
  48 + # Open reader
  49 + reader = csv.reader(open(filename, 'r'), delimiter=',')
  50 +
  51 + # Read rows
  52 + first = True
  53 + index = -1
  54 + for row in reader:
  55 +
  56 + # Get column index if first row
  57 + if first:
  58 + try:
  59 + index = row.index(parname)
  60 + except:
  61 + print 'ERROR: Parameter "'+parname+'" not found in list:'
  62 + for p in row:
  63 + print ' "'+p+'"'
  64 + raise NameError(parname)
  65 + #print index
  66 +
  67 + # Handle data rows
  68 + else:
  69 + values.append(float(row[index]))
  70 +
  71 + # Flag that first row has been passed
  72 + first = False
  73 +
  74 + # Create numpy array
  75 + a = np.array(values)
  76 + #print a
  77 +
  78 + # Return array
  79 + return a
  80 +
  81 +
  82 +# ======================== #
  83 +# Main routine entry point #
  84 +# ======================== #
  85 +if __name__ == '__main__':
  86 + """
  87 + Display pull histogram generated by cspull using matplotlib.
  88 + """
  89 + # Print usage information
  90 + usage = "Usage: pull_histogram filename parname [bins]"
  91 + if len(sys.argv) < 3:
  92 + print usage
  93 + sys.exit()
  94 +
  95 + # Extract parameters
  96 + filename = sys.argv[1]
  97 + parname = sys.argv[2]
  98 + if len(sys.argv) > 3:
  99 + nbins = int(sys.argv[3])
  100 + else:
  101 + nbins = 50
  102 +
  103 + # Read values from CSV file
  104 + values = read_pull(filename, parname)
  105 +
  106 + # Create histogram
  107 + n, bins, patches = plt.hist(values, nbins, range=[-4.0,4.0],
  108 + normed=True, facecolor='green')
  109 +
  110 + # Create expected distribution
  111 + y = mlab.normpdf(bins, 0.0, 1.0)
  112 + l = plt.plot(bins, y, 'r-', linewidth=2)
  113 +
  114 + # Set parname for plotting
  115 + name = parname
  116 +
  117 + # Set plot
  118 + plt.xlabel('Pull ('+name+')')
  119 + plt.ylabel('Arbitrary units')
  120 + plt.title(name)
  121 + plt.grid(True)
  122 +
  123 + # Show histogram
  124 + plt.show()
  125 +
0 126 \ No newline at end of file
... ...