2. Collection of useful tools¶

2.1. Overview¶

Below are a few programming examples to create tools that help in the analysis of experimental data such as:

Calculate a numerical derivative
Use the zoom tool of a plot to select a sub-set of data
Determine the position of a peak
How to fit a gaussian on a linear back ground.

Example data gauss_bkg.data. To use the examples below load the data as follows:

d = B.get_file('gauss_bkg.data')
xv = d['x']
yv = d['S']
dyv = d['dS']

2.2. Calculate a derivative numerically¶

When you have to calculate the propagation of measurement errors it is necessary to calculate partial derivatives. For some functions that can get quite complex and a simpler approach is to calculate this numerically. A simple approximation that is accurate to second order is the symmetric difference quotient for a function \(f(x)\) with \(h\) as step size:

(2.1)¶\[\frac{f(x + h) - f(x - h)}{2 h}\]

This can be programmed as follows:

# calculate the derivative of function f(x)
def dfdx(f, x, eps = 0.01):
   """
   return an estimate of the derivative of f at x.

   eps: determines the step size as a fraction of x
   """
   h = eps*x
   return (f(x + h) - f(x - h))/(2*h)

If you need a different step size use the keyword eps to change it by calling this function like dfdx(my_func, x, eps = 0.005), where my_func represents your function.

2.3. Select a sub-set of data interactively¶

Frequently you are confronted with the task of working with a subset of your data. This happens frequently when you have to determine the location or the height of a peak in your data. One possibility to do this is to plot all your data and then try to read the limits off your x-axis. This methods can become time-consuming an can lead to errors. Instead you can use the following method:

# get a selection using a plot
def get_view(x):
   """
   return a selection for data currently plotted

   x array of data that were plotted
   """
   x1,x2 = B.pl.xlim()
   s = (x1 <= x)&(x<=x2)
   return s

This can be used as follows: assuming xv, yv and dyv are your experimental data that you plotted using:

B.plot_exp(xv, yv, dyv)

By selecting the zoom tool of the plot zoom into a part of your plot. In the console window do now the following:

s = get_view(xv)
print(s)

You will see that s is an array containing True and False values. Un-zoom your plot by clicking the Home icon and plot the selected data:

B.plot_exp(xv[s], yv[s], dyv[s])

and you will see that only the data from the previous zoomed view that you selected using get_view are plotted in a different color. Selecting data in this way is very handy when you have to fit peaks and work on sub-sets of your data.

2.4. Fitting a Peak in a Histogram¶

Here is an example of a function that makes determining a peak position in a histogram easier. It is based on the same idea that you zoom into the peak you wan to fit and extract the plot limits which you then use in your fit.

First let’s load a Cobalt 60 spectrum (you can use co_60_cal.SPE as an example):

# for this load an example spectrum
sp = B.get_spectrum('co_60_cal.SPE')

The function below will fit a gaussian to the data shown and return the position in its error:

def peakpos(h):
   """
   h : histogram you would like to fit
   fit that part of the histogram shown
   you need to have executed h.plot() and zoomed in to
   the peak region you want to fit
   """
   # get current limits
   x1, x2 = B.pl.xlim()
   # fit a gaussian to the histgram between x1 and x2
   h.fit(x1,x2)
   # plot the fit
   h.plot_fit()
   # return the peak position value
   return h.mean.value, h.mean.err

You can use this as follows:

sp.plot()
# zoom into one of the peaks
p, dp = peakpos(sp)
# p and dp contains the fitted peak position and the fit it is plotted as well.

2.5. Peak Position by Fitting a Parabola¶

The idea behing this method is to approximate the top part of a peak by a parabola. The parameters of the parabola are determined by a least square fit. It is easy to determine the location of the extremum (maximum or minimum) of a parabola which would correspond to the position of your peak. This of course only works if you select only the top part of your observed peak. An example function is given below:

# determine the peak position by fitting a parabola to data and determine its
# extremum location

def get_peak(x,y,dy = None, plot = False):
   """
   Parameters
   ----------
   x : independent data
   y : exp.  data
   dy : error in exp. data

   Returns
   -------
   peak position, uncertainty in peak position

   """
   # fit parabola to data
   if dy is None:
       p = B.polyfit(x, y, order =2)
   else:
       p = B.polyfit(x, y, dy, order =2)
   #
   # y = a0 + a1*x + a2*x**2
   #
   a0 = p.par[0]
   a1 = p.par[1]
   a2 = p.par[2]
   # covariance matrix
   C = p.cov
   # parabola max/min
   xm = -a1/(2.*a2)
   # for errors calculate partial derivatives wrsp a1, a2
   dxm_da1 = -1./(2.*a2)
   dxm_da2 = a1/(2*a2**2)
   # calculate total error using the covariance matrix
   dxm2 = (C[1,1]*dxm_da1**2 + 2.*C[1,2]*dxm_da1*dxm_da2 + C[2,2]*dxm_da2**2)
   dxm = np.sqrt(dxm2)
   # if selected plot the fitted curve
   if (plot):
       B.plot_line(p.xpl, p.ypl)
   return (xm, dxm)

This function can be used as follows. Plot the data and zoom into your peak. Make sure that the left and right edges of your plot define the edges of the top part:

# zoom to the top part of the peak
# get the selection
s = get_view(xv)
# get the peak position and ask for a plot of the fit
xmax, dxmax = get_peak(xv[s], yv[s], dyv[s], plot = True)
# xmax : peak position
# dxmax : error in peak position

2.6. Fitting a Gaussian¶

Gaussians are frequently the standard function used to represent a peak. Below is an example how to you the genfit function to fit a gaussian to data. First you need to define the fit parameters. I used A for the peak height, x0 for the peak position and sigma for the width. The parameters are defined as follows:

# define fit parameters and initialize them

A = B.Parameter(1.,'A')  # initial value 1
x0 = B.Parameter(1.,'x0') # initial value 1
sigma = B.Parameter(1.,'sigma') # initial value 1

After that you define the fit function using these parameters:

# define gaussian function using parameters
def gauss(x):
  y = A()*np.exp(-((x-x0())/(2.*sigma()))**2)
  return y

Now you are ready for the fit-function:

# this is the main function
def gauss_fit(x, y, dy):
  """

  Parameters
  ----------
  x : array of indep. var. values
  y : array of dependent var. values
  dy : array of y errors

  Returns
  -------
  fit results

  """
  g = B.genfit(gauss, [A, x0, sigma], x, y, y_err = dy)
  return g

You can use this function now as follows:

fit = gauss_fit(xv[s], yv[s], dyv[s])
# plot the fitted line
B.plot_line(fit.xpl, fit.ypl)

A more compact and modern way of programming the same task would be to use a class:

# fit a gaussian by creating a class (object oriented programming example)

class gauss_fit:
    # initial values can be given when a fit object is created
    def __init__(self, A=1., x0 = 1., sigma = 0.1):
        # initialize the class and store initial values
        self.A = B.Parameter(A, 'A')
        self.x0 = B.Parameter(x0, 'x0')
        self.sigma = B.Parameter(sigma, 'sigma')
        self.fit_list = [self.A, self.x0, self.sigma]

    def gauss(self, x):
        y = self.A()*np.exp(-((x-self.x0())/(2.*self.sigma()))**2)
        return y

    def fit(self, x, y, dy):
        self.x = x
        self.y = y
        if dy is None:
            self.dy = np.ones_like(y)
        else:
            self.dy = dy
        self.fit = B.genfit(self.gauss, self.fit_list, self.x, self.y, y_err = self.dy)
    def plot_fit(self):
        B.plot_line(self.fit.xpl, self.fit.ypl)
# end of class definition
# the peak position in this case is the parameter x0()

This can be used as follows:

# create the fit object:
G = gauss_fit()
# set initial parameters
G.A.set(50.);G.x0.set(5.);G.sigma.set(.5)
# do the fit:
G.fit(xv[s], yv[s], dyv[s])
# plot the fit
G.plot_fit()

The last example shows how to fit a gaussian on a linear background. This is also programmed as a class:

# A gaussian on a linear back ground
class gauss_bkg_fit:
    def __init__(self, A=1., x0 = 1., sigma = 0.1, a0 = 0., a1 = 0.):
        # peak parameters
        self.A = B.Parameter(A, 'A')
        self.x0 = B.Parameter(x0, 'x0')
        self.sigma = B.Parameter(sigma, 'sigma')
        # back ground parameters
        self.a0 = B.Parameter(a0, 'a0')
        self.a1 = B.Parameter(a1, 'a1')
        # fit list, which parameters are to be varied
        self.fit_list = [self.A, self.x0, self.sigma,self.a0, self.a1]

    def gauss(self, x):
        # gaussian
        y = self.A()*np.exp(-((x-self.x0())/(2.*self.sigma()))**2)
        return y
    def bkg(self, x):
        # lin. back ground
        y = self.a0() + self.a1()*x
        return y
    def signal(self,x):
        # combine the 2 functions
        return self.gauss(x) + self.bkg(x)

    def fit(self, x, y, dy):
        self.x = x
        self.y = y
        if dy is None:
            self.dy = np.ones_like(y)
        else:
            self.dy = dy
        self.fit = B.genfit(self.signal, self.fit_list, self.x, self.y, y_err = self.dy)
    def plot_fit(self):
        B.plot_line(self.fit.xpl, self.fit.ypl)
# end of class definition

And this is used similarly to the previous example:

Gb = gauss_bkg_fit()
# set initial parameters
Gb.A.set(50.);Gb.x0.set(5.);Gb.sigma.set(.5)
# do the fit:
Gb.fit(xv, yv, dyv)
# plot the fit
Gb.plot_fit()

You can control which parameters are fit by changing the content of Gb.fit_list above. In all the above examples of fitting a gaussian the peak position is the value of the parameter X0. You can access the values of parameters as follows:

print(x0)
xp = x0.value
dxp = x0.err
name = x0.name

Table of Contents

Related Topics