Michael Mayer / ctools

Commit 21078497a6b15e097c0d3d16f9bdab90c4388d46

Authored by Jürgen Knödlseder
1 parent 396a007c

Add show_pull_evolution.py script. 

The show_pull_evolution.py shows the evolution of the pull distribution mean and standard deviation as function of the number of samples. It is mainly used to assess from how many samples on the pull distributions becomes relatively stable.
Inline Side-by-side
Showing 1 changed file with 124 additions and 0 deletions
examples/show_pull_evolution.py 0 → 100755
  View file @2107849
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script displays the evolution of the pull distributions as function
  4 +# of the number of samples. The pull distributions are generated using the
  5 +# cspull script. The pull is defined by
  6 +#
  7 +# (fitted value - simulated value) / estimated error
  8 +#
  9 +# and the pull histogram allows to assess possible biases in the fit
  10 +# parameters and errors.
  11 +#
  12 +# Required 3rd party modules:
  13 +# - matplotlib
  14 +# - numpy
  15 +#
  16 +# Copyright (C) 2015 Juergen Knoedlseder
  17 +#
  18 +# This program is free software: you can redistribute it and/or modify
  19 +# it under the terms of the GNU General Public License as published by
  20 +# the Free Software Foundation, either version 3 of the License, or
  21 +# (at your option) any later version.
  22 +#
  23 +# This program is distributed in the hope that it will be useful,
  24 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  25 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  26 +# GNU General Public License for more details.
  27 +#
  28 +# You should have received a copy of the GNU General Public License
  29 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  30 +#
  31 +# ==========================================================================
  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 + sys.stdout.write('ERROR: Parameter "'+parname+'" not found in list:\n')
  62 + for p in row:
  63 + sys.stdout.write(' "'+p+'"\n')
  64 + raise NameError(parname)
  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 +
  76 + # Return array
  77 + return a
  78 +
  79 +
  80 +# ======================== #
  81 +# Main routine entry point #
  82 +# ======================== #
  83 +if __name__ == '__main__':
  84 + """
  85 + Display pull evolution generated by cspull using matplotlib.
  86 + """
  87 + # Print usage information
  88 + usage = "Usage: show_pull_evolution filename parname"
  89 + if len(sys.argv) < 3:
  90 + sys.stdout.write(usage+"\n")
  91 + sys.exit()
  92 +
  93 + # Extract parameters
  94 + filename = sys.argv[1]
  95 + parname = sys.argv[2]
  96 +
  97 + # Read values from CSV file
  98 + values = read_pull(filename, parname)
  99 +
  100 + # Compute mean and rms
  101 + nsamples = len(values)
  102 + samples = []
  103 + mean = []
  104 + rms = []
  105 + for i in range(1,nsamples):
  106 + samples.append(float(i+1))
  107 + mean.append(np.mean(values[0:i]))
  108 + rms.append(np.std(values[0:i]))
  109 +
  110 + # Plot mean and rms evolution
  111 + plt.semilogx(samples, mean, 'r-', label="Mean")
  112 + plt.semilogx(samples, rms, 'b-', label="Std. deviation")
  113 +
  114 + # Set plot
  115 + plt.xlabel("Samples")
  116 + plt.ylabel("Value")
  117 + plt.title(parname)
  118 + plt.grid(True)
  119 +
  120 + # Set legend
  121 + plt.legend(loc="upper left")
  122 +
  123 + # Show histogram
  124 + plt.show()
... ...