Collection of useful tools =========================== 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 :download:`gauss_bkg.data <./example/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'] 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 :math:`f(x)` with :math:`h` as step size: .. math:: :label: der \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. 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. 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 :download:`co_60_cal.SPE <./example/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. 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 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 .. include:: include/links.rst