All these examples are in the examples subdirectory of the HippoDraw installation. To run the function_ntuple.py for example, just type the follow to your shell prompt.
> python -i function_ntuple.py
The -i option means use the name script as an initialization script for the Python session. After the script completes, one can type further commands in the Python shell. An alternative is to first start Python then type the following in the Python shell...
>>> execfile ( 'function_ntuple.py' ) >>>
To run all the examples, type ...
> python -i run_examples.py
All the examples start with a line to create an instance of the HippoDraw application like the one shown below.
from setPath import app, canvas
You can type this interactively as well. One can get on-line help for the HippoDraw extension module by using Python's built in help facility. Here's how ...
>>> import hippo >>> help (hippo)
Or for an individual component of the module like this...
>>> import hippo >>> help (hippo.Display)
The next section lists some typical tasks one might want to do from Python. For each task, an example is listed to illustrate how it is done.
NTuple creation by reading ASCII file:
NTuple creation by reading ROOT file: NTuple creation by reading FITS file: Finding NTuple names from a file: Adding rows to an NTuple: NTuple creation from Python lists: Create and use DataArray:The Python code is shown below.
"""
Demonstraton of static versus dynamic histograms. Also demonstrates
taking difference between two histograms and plotting it.
@author Paul F. Kunz <Paul_Kunz@slacs.stanford.edu>
"""
import random
from time import sleep
from setPath import app, canvas
from hippo import Display, NTuple, NTupleController
# Create static histogram and add it to canvas.
sthist = Display ( "Static Histogram" )
sthist.setTitle ( "Gaussian Distribution (static hist)" )
sthist.setLabel ( 'x', 'X' )
sthist.setRange ( 'x', 0, 100)
sthist.setBinWidth ( 'x', 1. )
canvas.addDisplay ( sthist )
# Create second static histogram and add it to canvas
sthists = Display ( "Static Histogram" )
sthists.setTitle ( "Gaussian Distribution (low statistics)" )
sthists.setLabel ( 'x', 'X' )
sthists.setRange ( 'x', 0, 100)
sthists.setBinWidth ( 'x', 1. )
canvas.addDisplay ( sthists )
# Create empty NTuple and set the column label.
# Setting the column labels sets the number of columns
ntuple = NTuple ( )
ntuple.setTitle ( 'Gaussian Distribution' )
ntuple.setLabels ( ('X value', ) )
# Register the NTuple so the Inspector can see it.
NTupleController.instance().registerNTuple ( ntuple )
# Create dynamic histogram and attach it to NTuple.
# It will automatically adjust its range, title, and labels
dyhist = Display ( "Histogram", ntuple, ('X value', ) )
canvas.addDisplay ( dyhist )
mean = 45
sigma = 10
gauss = random.gauss
# Generate some data, fill the static histograms and NTuple.
for i in range(10000):
x = gauss( mean, sigma )
sthist.addValues ( (x, ) )
ntuple.addRow ( (x, ) )
if i < 1000 :
sthists.addValues ( (x, ) ) # only fill with first 1000
# Print some statistics from static histogram
# Could do same for dynamic
datarep = sthist.getDataRep()
print "Histogram :"
print " Title : " + sthist.getTitle()
print " Entries : %i" % sthist.numberOfEntries()
print " Mean = %f" % datarep.getMean ( 'x' )
print " Rms = %f" % datarep.getRMS ( 'x' )
# Print the average X value on the display
canvas.selectDisplay ( sthist)
canvas.addTextRep ( sthist, 'averagex' )
canvas.selectDisplay ( dyhist)
canvas.addTextRep ( dyhist, 'averagex' )
# Get the contents of the bins as a DataArray
high = sthist.createDataArray()
low = sthists.createDataArray ()
# Take difference with high statistics one scaled down, and a column
# to the low one.
low[ 'diff' ] = high[ 'Entries / bin' ] / 10 - low[ 'Entries / bin' ]
# Create an XY plot to see the difference
xyplot = Display ( "XY Plot", low, ( 'X', 'diff', 'nil', 'Error' ) )
canvas.addDisplay ( xyplot )
low.register ('difference ntuple' )
print "Static histogram is on the left, dynamic one on the below it."
print "You can changing the binning of the dynamic one with the Inspector"
print "\nOn the right we have low statistics gaussian with difference"
print "between it and high statistics plotted below it"
The script first creates static histograms, sets their attributes, and adds it to the canvas. It then creates an empty NTuple, sets its title, and labels for the columns. Since there is only one string in the tuple of labels, the ntuple will have only one column. The dynamic histogram is then created, bound to the ntuple, and adding to the canvas.
The script generates some data. In the for loop, it fills one static histogram with all the generated data and the other one with only the first 1000 "events". It also adds a rows to the ntuple. Both histograms are initially empty and you can watch them grow when you run the script.
When the data generation is done, some statistical quantities are retrieved from the static histogram and printed. One can now use the Inspector Users Guide to change the attributes of the histograms. With the dynamic histogram one can change the size of the bins, but not with the static one. They are then added to overlay the plot.
Finally the we take the difference between the high and low statistics histograms. This is done by creating a DataArray which wraps an NTuple containing the contents of the bins. DataArray objects a column of data as a numarray using syntax of a Python dictionary. Numarray objects overload the arithmetic operators, so we can take the difference with vector algebra. Adding the created column to the small DataArray also has the syntax of adding an entry to a Python dictionary. An XY plot is then created using DataArray.
""" -*- mode: python -*-
Create an NTuple with power-law distribution and histogram it on log log scale
@author J. Chiang <jchiang@slac.stanford.edu>
"""
from setPath import app, canvas
import random, math
shoot = random.uniform
# Create empty NTuple of 4 columns that we can fill by row.
from hippo import NTuple, NTupleController
nt = NTuple ()
nt.setLabels ( [ 'x' ] )
xmin = 1.
xmax = 100.
gamma = 2.1
xpmax = math.pow ( xmax, 1. - gamma )
xpmin = math.pow ( xmin, 1. - gamma )
nsamp = 10000
for i in range(nsamp):
x = shoot(0, 1)
xx = math.pow ( x*( xpmax - xpmin ) + xpmin, 1./ (1. - gamma) )
nt.addRow ( (xx, ) )
from hippo import Display
hist = Display ( 'Histogram', nt, ( 'x', ) )
canvas.addDisplay ( hist )
# Set the x and y axis on log scale
hist.setLog ( 'x', True )
hist.setLog ( 'y', True )
# fit to power law function
from hippo import Function
datarep = hist.getDataRep ()
powerlaw = Function ( "PowerLaw", datarep )
powerlaw.addTo( hist )
powerlaw.fit()
After creating of the histogram, the display is set to be a log-log plot. When the X axis is set to log scale, the width of the bins are logarithmically increasing in size. They only appear uniform in width on the screen. The function used is a power law, which appears as straight line on this log-log plot.
""" -*- python -*-
This script demonstrates the creation of a displays with multiple data reps
Author: Paul_Kunz@slac.stanford.edu
$Id: datareps.py,v 1.4 2004/07/14 05:42:58 pfkeb Exp $
"""
from setPath import app, canvas
from hippo import NTupleController
# Create NTuple with its's controller.
ntc = NTupleController.instance()
nt = ntc.createNTuple ( 'aptuple.tnt' )
from hippo import Display
# Create a histogram
hist = Display ("Histogram", nt, ('Age', ) )
canvas.addDisplay( hist )
# Overlay another histogram.
hist.addDataRep ( 'Histogram', nt, ['Service'] )
reps = hist.getDataReps ()
reps[1].setLineStyle ( 'Dash' )
reps[1].setColor ( 'red' )
reps[0].setBinWidth ( 'x', 1. )
reps[1].setBinWidth ( 'x', 2. )
nt2 = ntc.createNTuple ( 'Z0Strip.data' )
# Create a strip chart display.
strip = Display ( "Strip Chart", nt2, ('Energy', 'Sigma' ))
canvas.addDisplay ( strip )
# Overlay strip chart with an XY plot.
strip.addDataRep ( "XY Plot", nt2, ['Energy', 'Sigma', 'nil', 'error'] )
# Get the second data representation to set its symbol and size
reps = strip.getDataReps ()
reps[1].setSymbol ( 'circle', 8. )
# Create 2D histogram as contour plot.
color = Display ( "Contour Plot", nt, ('Age', 'Cost' ))
color.setBinWidth ( 'x', 1.0 )
canvas.addDisplay ( color )
# Overlay it with 1D profile plot.
color.addDataRep ( 'Profile', nt, [ 'Age', 'Cost' ] )
print "Overlayed various kinds of data representations"
With HippoDraw one can overlay data representations of different kinds. One can even overlay two histograms with different bin widths as is shown in this example. One can also change the line style or symbol type as well as the color. How to do these is easily seen in the Python code.
""" -*- python -*-
This script adding functions and fitting. It also demonstrates
retreiving an ntuple from the histogram to do something with its
contents.
Author: Paul_Kunz@slac.stanford.edu
$Id: function_ntuple.py,v 1.3 2004/07/12 22:02:55 pfkeb Exp $nt
"""
from setPath import app, canvas
# Create NTuple with its controller so Inspector can see it.
from hippo import NTupleController
ntc = NTupleController.instance()
nt1 = ntc.createNTuple ( 'aptuple.tnt' )
from hippo import Display
hist = Display ( "Histogram", nt1, ("Cost", ) )
canvas.addDisplay( hist )
# Get the data representation so we can add function to it.
datarep1 = hist.getDataRep()
from hippo import Function
gauss = Function ( "Gaussian", datarep1 )
gauss.addTo ( hist )
# Get the function parameters and display them.
print "Before fitting"
parmnames = gauss.parmNames ( )
print parmnames
parms = gauss.parameters ( )
print parms
# Now do the fitting.
gauss.fit ( )
print "After fitting"
parms = gauss.parameters ( )
print parms
# Add another function.
gauss1 = Function ( "Gaussian", datarep1 )
gauss1.addTo ( hist )
# Do another fit, should fit to linear sum
gauss1.fit ()
# Add Chi-squared per d.f. display
canvas.addTextRep ( hist, 'Chi-squared' )
# Create an NTuple from the histogram.
# Calculate the residuals
result = hist.createNTuple ()
ntc.registerNTuple ( result )
coords = result.getColumn ( 'Cost' )
values = result.getColumn ( 'Density' )
res = []
for i in range ( result.rows ) :
x = coords[i]
diff = values[i] - gauss1.valueAt ( x )
res.append ( diff )
# Add a column and display it.
result.addColumn ( 'residuals', res )
resplot=Display ( "XY Plot", result, ( 'Cost', 'residuals', 'nil', 'Error' ) )
canvas.addDisplay ( resplot )
print "The histogram was fitted to the sum of two gaussians."
print "Then histogram bins were retrieved to calculate "
print "the residuals. These were then plotted as an XY Plot."
print 'One could have used the "Create residuals display" button on the'
print "Inspector, but that wouldn't have demonstrated anything."
After creating an NTuple from a ASCII file, a histogram is placed on the canvas. Then the data representation of it is retrieved, a function created and added to the representation. Next the function parameter names are retrieved along with the initial parameter values and printed. The parameter values are printed again after fitting. A second Gaussian function is added and the sum of the two is used for a fit. Then the chi-squared per degree of freedom is add to the plot.
The last part of the example illustrates retrieving information from the histogram. One asks the histogram to make an NTuple representation of itself. The NTuple is registered with the NTupleController so it will be visible from the Inspector. . One can thus access the histogram information by column. In the example, the residuals are calculated. The list of residuals are then added to the NTuple created by the histogram. Since this NTuple is of the same type as used to hold the initial data set, we can use it to make an XY plot as is shown in the example.
Everything done in this script could have been done with Inspector Users Guide. This includes creating of the residuals display as there is a button on the Function inspector to do this.
"""
Demonstrates making simple XY plot.
author Paul F. Kunz <Paul_Kunz@slac.stanford.edu>
"""
# Gets the HippoDraw application.
# import setPath
from setPath import app, canvas
from hippo import Display
# Create list of data
energy = [90.74, 91.06, 91.43, 91.50, 92.16, 92.22, 92.96, 89.24, 89.98, 90.35]
sigma = [ 29.0, 30.0, 28.40, 28.80, 21.95, 22.90, 13.50, 4.50, 10.80, 24.20]
errors = [ 5.9, 3.15, 3.0, 5.8, 7.9, 3.1, 4.6, 3.5, 4.6, 3.6]
# make a plot to test it.
xy = Display ( "XY Plot", [ energy, sigma, errors ],
['Energy', 'Sigma', 'nil', 'error' ] )
canvas.addDisplay ( xy )
xy.setTitle ( 'Mark II Z0 scan' )
print "An XY plot is now displayed. You can use the Inspector dialog"
print "to modify the appearance or fit a function to it."
After creating three Python lists of the data, one creates the display. The first argument to Display is the type of plot. The types available can be seen in the Data tabbed panel of the inspector. The second argument is a list of the data lists. This will have the effect of creating a NTuple and it will be visible from the Inspector. The third argument is also a Python list and it serves two purposes. It gives labels for the NTuple columns and it sets the binding options for the XY plot. The 3rd binding option is the X error, which is optional, so we use the "`nil'" string to indicate that binding is not to be used.
""" -*- python -*-
This script demonstrates various types of Displays.
Author: Paul F. Kunz <Paul_Kunz@slac.stanford.edu>
$Id: displays.py,v 1.5 2005/06/25 13:23:45 pfkeb Exp $
"""
from setPath import app, canvas
from hippo import NTupleController
# Read an ASCII NTuple file to create an NTuple.
ntc = NTupleController.instance()
nt1 = ntc.createNTuple ( 'aptuple.tnt' )
from hippo import Display
# Change the matrix from the default
canvas.setPlotMatrix ( 3, 4 )
# Create the displays.
hist = Display ("Histogram", nt1, ('Cost', ) )
canvas.addDisplay( hist )
profile = Display ( "Profile", nt1, ('Age', 'Cost' ))
canvas.addDisplay ( profile )
color = Display ( "Color Plot", nt1, ('Age', 'Cost' ))
color.setBinWidth ( 'x', 1.0 )
canvas.addDisplay ( color )
contour = Display ( "Contour Plot", nt1, ('Age', 'Cost' ))
contour.setBinWidth ( 'x', 1.0 )
canvas.addDisplay ( contour )
profile2D = Display ( "Profile 2D", nt1, ('Age', 'Cost', 'Service' ))
canvas.addDisplay ( profile2D )
profileC = Display ( "Profile Contour", nt1, ('Age', 'Cost', 'Service' ))
canvas.addDisplay ( profileC )
scatter = Display ( "Scatter Plot", nt1, ('Age', 'Cost' ))
canvas.addDisplay ( scatter )
nt2 = ntc.createNTuple ( 'Z0Strip.data' )
strip = Display ( "Strip Chart", nt2, ('Energy', 'Sigma' ))
canvas.addDisplay ( strip )
xy = Display ( "XY Plot", nt2, ('Energy', 'Sigma', 'nil', 'error' ))
canvas.addDisplay ( xy )
y = Display ( "Y Plot", nt2, ('Energy', ) )
canvas.addDisplay ( y )
xyz = Display ( "XYZ Plot", nt1, ('Age', 'Service', 'Cost' ))
canvas.addDisplay ( xyz )
print "Various type of displays were created"
The beginning of the script sets the form factor of how the plots will be displayed, name 3 across and 4 down. This was done so they would all fit on one page. Nothing else is new compared to the previous examples.
""" -*- mode:python -*-
Demo of reading ROOT file with function, cuts, and calculation
author: Paul F. Kunz <Paul_Kunz@slac.stanford.edu>
"""
from setPath import app, canvas
from hippo import RootController, Display
rc = RootController.instance()
filename = "surface_muons_4M_merit.root"
# Get the top level tree names from the file
ntuple_names = rc.getNTupleNames ( filename )
print "In this file, tree names are ", ntuple_names
# Create a data array from the first tree
ntuple = rc.createDataArray ( filename, ntuple_names[0] )
# Pprint some information about this tree
print "Number of columns = ", ntuple.columns
labels = ntuple.getLabels()
print "First ten olumn labels are ... ", labels[:10]
print "Number of rows = ", ntuple.rows
# Create a histogram for one of the columns and add it to the canvas
hist = Display ( "Histogram", ntuple, ('TkrEnergy', ) )
canvas.addDisplay ( hist )
# Set the Y axis on log scale of better viewing
hist.setLog ( 'y', True )
# Add a cut from data in another column
from hippo import Cut
hits_cut = Cut ( ntuple, ('TkrTotalHits',) )
canvas.addDisplay ( hits_cut )
hits_cut.setLog ( 'y', True )
hits_cut.addTarget ( hist )
hits_cut.setCutRange ( 4, 110, 'x' )
# Change the range of the displayed data
hist.setRange ( 'x', 40, 700 )
# fit a function to the histogram
from hippo import Function
datarep = hist.getDataRep ()
exp1 = Function ( "Exponential", datarep )
exp1.addTo ( hist )
exp1.fit ()
# Print the results of the fit
pnames = exp1.parmNames ()
print pnames
parms = exp1.parameters ()
print parms
# add another function to the histogram and fit the linear sum
exp2 = Function ( "Exponential", datarep )
exp2.addTo ( hist )
exp1.fit() # always fit to linear sum
# Build an array which is the sum of two columns and add it to the ntuple
label = "Raw sum"
ntuple [ label ] = ntuple [ 'TkrEnergy' ] + ntuple [ 'CalEnergySum' ]
# Histogram the new column
sum_hist = Display ( 'Histogram', ntuple, (label, ) )
canvas.addDisplay ( sum_hist )
sum_hist.setLog ( 'x', True )
sum_hist.setLog ( 'y', True )
# Compare it with the official sum of energies
merit_sum = Display ( 'Histogram', ntuple, ( 'EvtEnergySumOpt', ) )
canvas.addDisplay ( merit_sum )
merit_sum.setLog ( 'x', True )
merit_sum.setLog ( 'y', True )
sum_hist.setRange ( 'x', 1.0, 1e+06 )
One uses the single instance of the RootController object to read the file. First one gets the names of the top level TTree which HippoDraw interprets as NTuple names.
In the next step, one could have created an NTuple with the method createNTuple and the same arguments. It would be implemented by the RootNTuple class. This class only reads the data into memory as it is needed, so it is quite efficient even for files that have a large number of columns.
Instead of a NTuple, one creates a DataArray. It contains a NTuple and behaves like one for most purposes but gives us additional flexibility as will be seen near the end of the example.
The next few lines shows how to get the number of rows and columns as well as the labels for the columns. The script then creates a histogram and adds it to the canvas. Since the distribution falls off rather quickly, the Y axis is set to a log scale.
The script then applies a cut and adds the display of the cut variable to the canvas. The cut distribution also falls rather quickly to its Y axis is also put on a log scale.
The next step gets a handle on the data representation, i.e. the actual histogram, in order to fit a function to it. Not satisfied with the fit, a second function is applied. Note that when two functions are applied, a fit is done to the linear sum of the functions.
The last part of the script shows the power of using a DataArray. In one line the script adds two columns together, row by row, and adds the result to the DataArray. Let's take it step by step to understand it. The expression
ntuple [ 'TkrEnergy' ]
shows that the DataArray behaves like a Python dictionary in that we can access a column by a key. The key is the label of the column. It returns a numarray. The expression
ntuple [ 'TkrEnergy' ] + ntuple [ 'CalEnergySum' ]
thus accesses two columns and adds them together because the numarray class overload the operator + to create a new numarray whose values are the sum of the others, element by element.
Finally, a column is added with the expression
ntuple [ label ] =
where again the DataArray behaves like a Python dictionary. The right hand side of the assignment must be a numarray, which it is in this case. The column is not added to the ROOT file, it is only logically added in memory.
Treating columns of the data source as vectors and doing vector algebra on them can be considerably faster than writing for loops. The code is much easier to write and understand as well.
The instance of the FitsController is used to read the FITS file. FITS files can contain multiple images and tables. You can use the FitsController to find the names of the Header/Data Units (HDU) as shown. The first HDU is always an image, so the tabluar data is probably the second HDU. The FitsController creates FitsNTuple object given the filename and the name of the HDU.
The FitsController and FitsNTuple classes are written in C++ and use the standard CFITSIO package for reading the FITS file.
""" -*- python -*-
This script demonstrates using FITS file as data source.
Author: Paul F. Kunz <Paul_Kunz@slac.stanford.edu>
$Id: mainpage.py,v 1.3 2005/06/24 17:20:51 pfkeb Exp $
"""
from setPath import app, canvas
from hippo import FitsController
# Read an FITS table from file to create an NTuple.
ntc = FitsController.instance()
#
# Find out what the names of the HDUs in the file.
#
hdus = ntc.getNTupleNames ( 'SimEvents.fits' )
print "Names of the Header/Data Units are..."
print hdus
events = ntc.createNTuple ( 'SimEvents.fits', hdus[1] )
from hippo import Display
# Create the displays.
hist = Display ("Histogram", events, ('energy', ) )
hist.setLog ( 'x', True )
hist.setLog ( 'y', True )
canvas.addDisplay( hist )
color = Display ( "Color Plot", events, ('GLON', 'GLAT' ))
color.setRange ( 'z', 0.5, 150. )
color.setLog ( 'z', True )
canvas.addDisplay ( color )
contour = Display ( "Contour Plot", events, ('time', 'GLON' ))
canvas.addDisplay ( contour )
profile = Display ( "Profile", events, ('time', 'energy' ))
canvas.addDisplay ( profile )
print "Created the plots shown on the HippoDraw main page."
This example also illustrates how to directly manipulate the fitting process, using various components. In future releases, some of these components could also be written in Python.
Finally, this example uses a DataArray to hold the data set.
""" -*- mode:python -*-
Demonstrates implement function in Python derived from the C++
abstract base class. Also demonstrates use of the fitter from Python,
without using the GUI
@author Paul F. Kunz <paul_kunz@slac.stanford.edu>
"""
#
#$Id: fitting.py,v 1.4 2005/05/31 14:40:33 pfkeb Exp $
#
from setPath import app, canvas
import random, numarray
slope = 2.
intercept = 1.
sigma = 2.
#
# Define a class in Python that is derived from C++ base class
#
from hippo import FunctionBase
class Linear ( FunctionBase ) :
def __init__ ( self, other = None ) :
if other : # copy constructor
FunctionBase.__init__( self, other )
else : # default constructor
FunctionBase.__init__( self )
self.initialize ()
def initialize ( self ) :
self.setName ( 'linear' )
self.setParmNames ( [ 'Intercept', 'Slope' ] )
self.setParameters ( [ intercept, slope ] )
def valueAt ( self, x ) :
parms = self.getParameters ()
return x * parms[1] + parms[0]
def derivByParm ( self, i, x ) :
if i == 0 :
return 1.0
if i == 1 :
return x
gauss = random.gauss
# Generate some data using numarray
npts = 20
x = numarray.arange(npts)
y = numarray.array( [gauss(slope*xx + intercept, sigma) for xx in x] )
xerr = numarray.ones ( npts )
yerr = sigma*numarray.ones(npts)
# Create a DataArray and register it so that the Inspector can display it.
# The DataArray is will wrap a DataSource of type NTuple.
from hippo import DataArray
da = DataArray ( 'NTuple' ) # use NTuple to store data
da.register ( 'randomized line' ) # name the data source
# Fill the contents by adding named columns in the order expected by
# the fitter (more later)
da['x'] = x
da['y'] = y
da['xerr'] = xerr
da['yerr'] = yerr
# Create an XY plot, ignoring xerr, and add it to the canvas.
from hippo import Display
disp = Display ('XY Plot', da, ('x', 'y', 'nil', 'yerr' ) )
canvas.addDisplay ( disp )
# Get the name of the available fitters from the factory
from hippo import FitterFactory
factory = FitterFactory.instance()
fitters = factory.names()
print fitters
# Create a fitter, in this case Minuit Migrad algorithm with Chi Sq
migrad = factory.create ( fitters[1] )
print migrad.name()
# create a function
f = Linear ()
print f.name()
print f.parmNames()
# Get the FCN so we can set it up
fcn = migrad.getFCN ()
fcn.setFunction ( f )
fcn.setDataSource ( da.dataSource() )
fcn.setUseErrors ( True )
# Since we are using Minuit, we can set limits on the parameters. If
# we wanted to do that we would uncomment the following.
#migrad.setLimits ( 0, -2., -1. )
# Minimize the Chi Sq
migrad.minimize()
# Print the resulting parameters
print f.getParameters()
# Add the function to the display of the data
disp.addFunction ( f )
# Add the function to the factory so it can be used from the GUI
# One might just do this instead of using the minimizers directly
from hippo import FunctionFactory
ff = FunctionFactory.instance()
ff.add ( f )
# Create an DataArray with columns not in order expected by the fitter
da2 = DataArray ( 'NTuple' ) # use NTuple to store data
da2.register ( 'randomized line 2' ) # name the data source
# Fill the contents by adding named columns in the order expected by
# the fitter (more later)
da2['y'] = y
da2['x'] = x
da2['yerr'] = yerr
da2['xerr'] = xerr
# Teach the FCN the columns to use, ignoring xerr. The `1' indicates
# the dimension of the coordinate space
fcn.setDataSource ( da2.dataSource(), 1, [ 1, 0, -1, 2 ] )
f.setParameters ( [ intercept, slope ] )
migrad.minimize ()
print f.getParameters ()
This example also illustrates how to directly manipulate the fitting process, using features, some of which are only available with the Minuit fitter.
""" -*- mode:python -*-
This example shows how to do 2 dimension fitting by talking to the
fitter and its FCN directly. It also shows the optional setting
for the fitting process that are only available with a Minuit based
fitter.
author Paul F. Kunz <Paul_Kunz@slac.stanford.edu>
"""
#
#$Id: fitting2.py,v 1.12 2005/06/17 22:11:27 pfkeb Exp $
#
## from setPath import *
## import setPath
# The following is just for developers who need to test their built
# version before it gets installed.
build_path = '../../hippodraw-BUILD' # top level of build directory
module_path = build_path + '/python/.libs' # path to hippomodule.so
import sys
sys.path.insert ( 0, module_path )
# end only for developers
import random, numarray
import hippo
from hippo import FunctionBase
#
# Declare a function that take is 2 dimensional.
#
class Planar ( FunctionBase ) :
def __init__ ( self, other = None ) :
if other :
FunctionBase.__init__( self, other )
else :
FunctionBase.__init__( self )
self.initialize ()
def initialize ( self ) :
self.setName ( 'Planar' )
self.setParmNames ( [ 'Interscept', 'SlopeX', 'SlopeY' ] )
self.setParameters ( [ 0., 1., 1. ] ) # initial parameters
def valueAt ( self, x, y ) :
parms = self.getParameters ()
return parms[0] + parms[1] * x + parms[2] * y
#
# Generate some data.
#
slope = 2.
intercept = 1.
mean = 10.
sigma = 4.
gauss = random.gauss
def plane ( x, y ) :
return x + y
npts = 20
x = numarray.arange(npts*npts)
y = numarray.arange(npts*npts)
z = numarray.arange(npts*npts)
zerr = numarray.arange(npts*npts)
k = 0
for i in range(npts) :
for j in range(npts) :
x[k] = xx = i
y[k] = yy = j
z[k] = plane ( xx, yy ) + gauss ( 0., 2. )
zerr[k] = sigma
k += 1
#
# Store the data in NTuple
#
nt = hippo.NTuple () # empty ntuple
nt.addColumn ( 'x', x ) # add column
nt.addColumn ( 'y', y )
nt.addColumn ( 'z', z )
nt.addColumn ( 'zerr', zerr )
#
# Get the list of available fitters so one can verify that Minuit is available.
#
from hippo import FitterFactory
factory = FitterFactory.instance()
fitters = factory.names()
print fitters
#
# Create a Minuit Migrad fitter
#
migrad = factory.create ( fitters[1] )
print migrad.name()
#
# Create an instance of our custom function and verify its names
#
f = Planar ()
print f.name()
print f.parmNames()
#
# Get the objective function so we can set the function and data to be used
#
fcn = migrad.getFCN ()
fcn.setFunction ( f )
#
# Tell the FCN how to bind to the NTuple
# columns are x y z zerr
fcn.setDataSource ( nt, 2, [ 0, 1, 2, -1, -1, 3 ] ) # no xerr or yerr
fcn.setUseErrors ( False ) # do not use errors on z value
migrad.minimize()
chisq = fcn.objectiveValue ()
print "ChiSq =", chisq
fcn.setUseErrors ( True ) # use error on z value.
chisq = fcn.objectiveValue ()
print "ChiSq =", chisq
#
# get a copy of thie parameter values
#
parms = list(f.getParameters())
print parms
#
# Set a parameter value and have it fixed in the fit
#
parms[1] = 1.0
f.setParameters ( parms )
migrad.setFixedFlags ( ( 0, 1, 0 ) )
migrad.minimize ()
print f.getParameters()
#
# Try setting limits, only available with Minuit fitter
#
migrad.setLimits ( 2, 0.99, 1.01 )
migrad.minimize()
print f.getParameters ()
# We don't need to display the data
# If we do, we would add lines like the commented out one that follow
#
## from hippo import Display
## disp = Display ( "XYZ Plot", nt, ( 'x', 'y', 'z' ) )
## canvas.addDisplay ( disp )
1.4.3