# Python Workshop I: Fitting a single symmetric peak

Time to play with real data! In this workshop, we will fit a peak with a single symmetric function. The script for the entire analysis is given below. Copy the whole thing into an editor and save it as ‘name_of_file.py’. I have used the Raman spectrum of silicon centered at the dominant peak (520 cm^{-1}). This data is available here. Please download it and put it into the same directory as the script. The code is pretty self-explanatory but here is an outline of the analysis routine:

**1. Import data**

**2. Background removal.** In this section, the regions outside the peak are isolated and fit to a linear background. To do this, we first describe the regions of the required data. The variables* ‘ind_bg_low’* and *‘ind_bg_high’* are actually boolean arrays created by the specified limits. So, *‘ind_bg_low’* will be ‘True’ for all *x* greater than the minimum value of *x *and less than *450*. These two variables are then used as ‘indices’ to splice the full *x* and *y* datasets.

**3. Fitting.** As with any optimization routine, one needs to define a function that can be fitted to the data. In this example, we have chosen a Lorentzian peak with three parameters that can be fitted: the half-width at half-maximum (hwhm), peak center and the peak intensity. These parameters are assigned to a list of parameters called *‘p’*. The items of the list, e.g. item 1(index=0), correspond to a given parameter. We have defined the list to have the parameters in this order:* hwhm, peak center, intensity*. In addition to defining a function, we also need to define a minimisation criterion. The* leastsq* function minimises the difference between the y values of the data and the fit. The function* ‘residuals’ *is defined such that it returns this difference and hence, it is the function that is given to the *leastsq* function call. The first part of the results outputted by the *leastsq* function contain the parameters of best fit. These parameters are then put into the variable *‘best_parameters’*. The fit is evaluated by substituting *‘best_parameters’ *into the Lorentzian function.

**4. Plotting data and fit**

This particular routine is an introduction to simple scenarios of each step. In the next workshop, we will look at more sophisticated background correction, fitting and plotting options.

######################################################################### ####################### IMPORTING REQUIRED MODULES ###################### import numpy import pylab from scipy.optimize import leastsq # Levenberg-Marquadt Algorithm # ######################################################################### ############################# LOADING DATA ############################## a = numpy.loadtxt('silicon_Raman.txt') x = a[:,0] y = a[:,1] ######################################################################### ########################### DEFINING FUNCTIONS ########################## def lorentzian(x,p): numerator = (p[0]**2 ) denominator = ( x - (p[1]) )**2 + p[0]**2 y = p[2]*(numerator/denominator) return y def residuals(p,y,x): err = y - lorentzian(x,p) return err ######################################################################### ######################## BACKGROUND SUBTRACTION ######################### # defining the 'background' part of the spectrum # ind_bg_low = (x > min(x)) & (x < 450.0) ind_bg_high = (x > 590.0) & (x < max(x)) x_bg = numpy.concatenate((x[ind_bg_low],x[ind_bg_high])) y_bg = numpy.concatenate((y[ind_bg_low],y[ind_bg_high])) #pylab.plot(x_bg,y_bg) # fitting the background to a line # m, c = numpy.polyfit(x_bg, y_bg, 1) # removing fitted background # background = m*x + c y_bg_corr = y - background #pylab.plot(x,y_bg_corr) ######################################################################### ############################# FITTING DATA ## ########################### # initial values # p = [5.0,520.0,12e3] # [hwhm, peak center, intensity] # # optimization # pbest = leastsq(residuals,p,args=(y_bg_corr,x),full_output=1) best_parameters = pbest[0] # fit to data # fit = lorentzian(x,best_parameters) ######################################################################### ############################## PLOTTING ################################# pylab.plot(x,y_bg_corr,'wo') pylab.plot(x,fit,'r-',lw=2) pylab.xlabel(r'$\omega$ (cm$^{-1}$)', fontsize=18) pylab.ylabel('Intensity (a.u.)', fontsize=18) pylab.show()

*Post contributed by Shrividya Ravi, VUW*

Thanks for your explanation. But what if I don’t have a numerical set of data but a microscope image with a laser groove in it?

First I did otsu-segmentation and output the intensity profile of the groove with .

profile = im0.sum(axis=0) #

plt.figure()

plt.plot(profile)

I then get an image with a peak. How to determine its FWHM?

Thanks, Sebastian

Hi Sebastian,

Please post here if you could find how to deteremine FWHM.

@James, very useful documentation. Thanks.

-Mithil

Hi Sebastian,

That’s a good question! I don’t know much about this, however there is a post on Stack Overflow which may help you. Pylab has some built in image processing functions, so it’s definitely possible with a little work.

http://stackoverflow.com/questions/14154233/image-analysis-curve-fitting

Hope that helps.

Thanks for this great example!

Just a (maybe not so minor) heads-up: The numerator of the Lorentzian in the example above is squared, but according to both MathWorld and Wikipedia it shouldn’t be.