Subsections

3.31 Fitting data. The HFit and HFitter classes

3.31.1 Preparing a fit

Let us first initialize the curve fitting package by creating an instance of the HFitter class. Then we will print available methods:

>>> from jhplot  import *
>>> f=HFitter()
>>> print f.getFitMethod()
['uml', 'leastsquares', 'cleverchi2', 'chi2', 'bml']

The most popular method is chi2 (this is the method which is used by the default when you create a HFitter() objects without passing any argument)

The package contains several pre-build functions which can be used for fitting. You can print the function catalogue as:

>>> from jhplot  import *
>>> f=HFitter()
>>> print f.getFuncCatalog()
['e', 'g', 'g2', 'lorentzian', 'moyal', \
 'p0', 'p1', 'p2', 'p3', 'p4', 'p5', \
 'p6', 'p7', 'p8', 'p9', 'landau', 'pow']

Table :autorefBuild-in fit functions of the HFitter class.3.1 shows the definitions.


Table: Build-in fit functions of the HFitter class.

Name
Definitions

e
Exponential
g Gaussian
g2 Double Gaussian
lorenzian Lorenzian
moyal Moyal
landau Landau
pow power-law, $a*(b-x)^c$
P0 $y=a$
P1 $y=a+b*x$
P2 $y=a+b*x+c*x^2$
Pn polynomial of n$^{th}$ order

Let us create a Gaussian function using the predefined key 'g' and print all attributes of this function:

>>> f.setFunc("g")
>>> func=f.getFunc()
>>> print func
BaseModelFunction:  Title=g, name=g
Dimension: 1, number of parameters: 3
Codelet String: codelet:g:catalog
Variable Names: x0,
Parameters: amplitude=1.0, mean=0.0, sigma=1.0,
Provides Gradient: true,  Provides Parameter Gradient: true,
Provides Normalization: false

You can also set parameter to a specific value. For example, one can set the mean of a Gaussian to 10 using this method:

>>> func.setParameter("mean",10)

One can also set all three parameters by using the method above:

>>> func.setParameters([100,10,1])

which sets the amplitude, mean and sigma to 100,10,1. One can build a Gaussian function plus a second-order polynomial using the string "g+p2", i.e.:

>>> f.setFunc("g+p2")
>>> func=f.getFunc()
>>> print func.title()+' has: ', func.parameterNames()
'amplitude', 'mean', 'sigma', 'p0', 'p1', 'p2'

One can also build a function from a script:

Let us illustrate how to build a parabolic function on-fly:

>>> f.setFunc('parabola',1, 'a*x[0]*x[0]+b*x[0]+c','a,b,c')
>>> func=f.getFunc()
>>> print func.title()+' has: ', func.parameterNames()
parabola has:  ['a', 'b', 'c']
>>> print func.codeletString()
codelet:parabola:verbatim:jel : 1 : a*x[0]*x[0]+b*x[0]+c : a,b,c :

The setFunc() method takes several arguments: function title, number of dimensions, a string representing the function, and the names of the free parameters. The variable of this function is "x[0]". In case, say 2D functions, you can add second variable as "x[1]".

One can set the parameter value as:

>>> f.setPar('mean',10)
This sets the mean to 10, assuming that setFunc('g') was run first as shown previously.

During the minimisation, one can restrict parameter variations to a certain range. This can be done with the method:

>>> f.setParRange('mean', -10, 10)

This sets the range [-10,10] of variations of the Gaussian means.

Next, we will try to fit our data with the prepared function. This part is easy: if you have P1D, H1D, H2D or PNDm just execute the method fit(data), where the data is one of the objects mentioned above:

>>> f.fit(data)

In case of 1D data, you may also restrict the fitting range with the method setRange(min,max). After the fit, you can get the resulting function and other characteristics as. The resulting function can be obtained as:

>>> func=f.getFittedFunc()

Another way to retrieve the results is to use this method:

>>> result=f.getResult()

So, what the object result can do? Well, this is a rather sophisticated class which keeps essentially everything you need to retrieve the results. Below we will list the most important methods and comment some of them:

}
>>> result.fittedParameters()
>>> result.errors()
>>> result.constraints()
>>> result.covMatrixElement(0,1)
>>> result.engineName()
>>> result.errors()       # parabolic errors
>>> result.errorsMinus()  # error -
>>> result.errorsPlus()   # error -
>>> result.fitMethodName()
>>> result.fitStatus() # status
>>> result.fittedFunction()
>>> result.fittedParameter("mean")
>>> result.fittedParameterNames()
>>> result.fittedParameters()
>>> result.ndf()          # number of degrees of freedom
>>> result.quality()      # chi2/ndf for Chi2 method

Below we give a complete example:

from jhplot  import *
from java.util import Random

f=HFitter()
print f.getFuncCatalog()
f.setFunc('g')
f.setPar('amplitude',50)

h1 = H1D('Data',50, -4, 4)
h1.setPenWidthErr(2)
h1.setStyle("p")
h1.setSymbol(4)
h1.setDrawLine(0)

r = Random()
for i in range(1000):
      h1.fill(r.nextGaussian())

c1 = HPlot('Canvas')
c1.visible()
c1.setAutoRange()
c1.draw(h1)

f.setRange(-4,4)
f.fit(h1)
ff=f.getFittedFunc()

# get fitted results
r=f.getResult()
Pars   = r.fittedParameters()
Errors = r.errors()
Names  = r.fittedParameterNames()
print "Fit results:"
for i in range(ff.numberOfParameters()):
  print Names[i]+" : "+str(Pars[i])+" +- "+str(Errors[i])


mess='χ^{2}/ndf='+str(round(r.quality()*r.ndf()))
mess=mess+" / "+str(r.ndf())
lab=HLabel(mess, 0.12, 0.69, 'NDC')
c1.add(lab)

# draw function
f2 = F1D("Gaussian",ff,-4,4)
f2.setPenWidth(3)
c1.draw(f2)
print 'Quality=',r.quality(), ' NDF=',r.ndf()

Check also examples : hfitter1.py and hfitter2.py.

3.31.2 Fitting data using FreeHEP libraries

Look at an example in histo_fit.py. This script fills a histogram, creates factories and does a fit with a Gaussian function.

3.31.3 HFit class

One can fit data (H1D histograms or P1D data holders) using the Fit Panel based on the HFit class. One can look at the example histo_fit1.py.

>>># .. define c1 as HPLot object and h1 as a H1D histogram and fill it 
>>>a=HFit(c1,"c1",h,"h1")  # start HFit panel
After execution of HFit(c1,"c1",h,"h1") (the arguments should be adjusted to reflect the variable names for HPlot and H1D objects), a pop-up fit panel will be shown. Note that we pass not only the objects, but also the variable names for the HPlot and H1D class. One can also pass a P1D object instead of a histogram H1D. One should pass the variable names if one needs to generate JAIDA source code automatically in future.

The fit panel allows the user to set up the predefined fit functions, fit method and perform the fit. The user can also set up the initial parameters for the fit function (use the Settings button). Once the fit is acceptable, the used can generate the output fit parameters, which will be inserted directly to the jeHEP editor in form of a Jython code. The user can also automatically generate the JAIDA source code which corresponds to the fit. Once inserted, it can further be corrected using the jeHEP editor. Then HFit() statement can be removed and the user can run the generated JAIDA code manually.

One can also specify a user-defined 1D function and add it to the HFit Panel. This is an example:

>>># .. define c1 as HPLot object and h1 as a H1D histogram and fill it
>>>a=HFit(c1,"c1",h,"h1")  # start HFit panel
>>>a.addFunc("User1", "Tool tip", "a*x[0]+b","a,b") # add "User1" function
After executing this scrip, a new function User1 will be added to the fit panel.