I have data with peaks on some background, for example:

enter image description here

The two prominent peaks at ~390 and ~450, as well as the much smaller peak at ~840. What are some options to programmatically find the position (i. e. the x-coordinate) of such peaks using Python/SciPy?

  • $\begingroup$ what you're asking for is pretty ad hoc and in the eye of the beholder. So maybe calculate the forward and backwards derivative at each point, if they have different sign then you have a peak and can store the value and location. But to winnow the small peaks (that you want to ignore) from the larger ones, maybe take a long windowed average (or function fit) and only keep the peaks more than 50% away from the average/function. You could play with the percentage till you get roughly what you want. At the end you'll get the value and position of every peak more than x% from the windowed average. $\endgroup$
    – EMP
    May 28 '19 at 21:52
  • 2
    $\begingroup$ docs.scipy.org/doc/scipy/reference/generated/… $\endgroup$
    – GoHokies
    May 29 '19 at 11:34
  • $\begingroup$ Just out of curiosity, what sort of data is this? NMR? Photo-emission? $\endgroup$ May 29 '19 at 15:01
  • $\begingroup$ @mathewgunther It's the energy of gamma rays from radioactive decay (some cobalt isotope). x-axis is the energy and y-axis is the number of events with that energy. $\endgroup$
    – 0x539
    May 30 '19 at 11:29
  • $\begingroup$ Do the peaks of interest have a closed form analytic expression? (e.g. Gaussian, Lorenzian, Voigt, etc...) Also, does the background and broad peaks below 380 have a closed form analytic expression? Finally, what sort of x-axis peak position resolution is needed? $\endgroup$ May 30 '19 at 18:19

Fitting the peaks of gamma spectra is a typical task in non-destructive analysis of spent fuel or neutron activation analysis. Since these applications are already "quite old", there is some standard software available, like Genie 2000. A paper Evaluation of Peak-Fitting Software for Gamma Spectrum Analysis from 2015 compares a number of these tools. However, most of them are commercial.

All methods are based on the concept of "base functions", typically Gaussian functions, that are used in a linear combination. The two main problems are then to estimate the number of peaks and to estimate their locations. Both are done using prior knowledge: what does the user think is present in the sample? Once this is done, a non-linear optimisation algorithm is used (Levenberg-Marquardt) to minimize the error between the fit and the measured data. But, as mentioned by EMP in a comment, a lot of this is in the eye of the beholder. Most commercial software have a GUI that allows the user to interact with the code to "tinkle until convergence" (where convergence means "it fits according to my human eye and according to what radioactive species I think I have in my material").

A paper describing a non-linear fitting procedure using the Levenberg-Marquardt method can be found here. This will give you some of the mathematical background. This paper "Peak fitting and identification software library for high resolution gamma-ray spectra" also gives some mathematical background. An older report describing this problem can be found here. Based on the latter and the available solvers in SciPy, you should be able to roll your own code.


That looks like a spetrum (of x-rays (?)). If you have any a-priori information about where you expect your peaks to land, you might give your fitting algorithm that information. E.g.: "fit a gauß peak between 390 and 415". That's a little ugly because you'd have to hand-code those ranges, but at least you can re-use it on scans later on automatically.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.