DMelt:DataAnalysis/9 Peak Finders

Table of contents

Peak finders

DataMelt has several Java algorithms for identifications of peaks in data. Let us consider a peak identification algorithms in arrays of data. It is important for peak identification in time series, mass spectra and other areas.

Given a spectrum and search parameters, performs a digital filter peak search as specified in V. Hnatowicz et al in Comp Phys Comm 60 (1990) 111-125. Setting the sensitivity to a typical value of 3 gives a 3% chance for any peak found to be false.Maximum separation in sigma between peaks is 1.3.

"""
Peak finder algorithm. Given a spectrum and search parameters, performs a digital filter peak search as specified in V. Hnatowicz et al in Comp Phys Comm 60 (1990) 111-125. Setting the sensitivity to a typical value of 3 gives a 3% chance for any peak found to be false.Maximum separation in sigma between peaks is 1.3.
See: https://datamelt.org/api/doc.php/jhpro/tseries/PeakFinder
"""
 
from jhplot import *
from jhpro.tseries import *

# create input data
spectrum=[10,10,10,10,10,10,10,40,10,10,10,10,10,50,10,11,10,10,10]
p=PeakFinder("test", spectrum, 1.2, 1.2)
mu=p.getPeaks()
print mu.tolist()

You can also find peak finder based on the algorithms described in the paper "A non-parametric peak finder algorithm and its application in searches for new physics" by S.Chekanov, M.Erickson, arxiv.org:1110.3772 "Advances in High Energy Physics", vol. 2013, Article ID 162986. This peak identification algorithm is implemented in Python:

"""
Peak finder algorithm based on
A non-parametric peak finder algorithm and its application in searches for new physics
S.Chekanov, M.Erickson
E-print: arxiv.org:1110.3772 ANL-HEP-PR-11-63
"Advances in High Energy Physics", vol. 2013, Article ID 162986, 4 pages, 2013. doi:10.1155/2013/162986.
http://www.hindawi.com/journals/ahep/contents/
http://www.hindawi.com/journals/ahep/2013/162986/ 
"""
 
from jhplot import *
from java.awt import Color
import sys
sys.path.append(SystemDir+'/python/packages/npfinder/')
from npfinder import GetData, FindPeaks, MakeGraph

c1 = HPlot("Canvas",600,400)
c1.visible()
c1.setLegend(0)
c1.setGTitle("Automatic peak identification")
c1.setRange(50,1000,1,1000)
# c1.setAutoRange()
# c1.setLogScale(1,1)

h=H1D("data",50,0,1000)
h.fillGauss(30000,0,300) # background 
h1=h.copy()
h.fillGauss(1000,400,20) # peak 1
h.fillGauss(500,700,10)  # peak 2 
data = GetData(h)
# 3-sigma approach
minSig=3
peaks = FindPeaks(data, 0.6, minSig) 
print 'RESULTS OF PEAK SEARCH:\n'
for i in range(len(peaks)):
    print 'Peak Number:', peaks[i].GetPeakNumber(), 'Peak Start:', peaks[i].GetPeakStart().x, 'Peak End:', peaks[i].GetPeakEnd().x, 'Peak Stat Sig:', peaks[i].GetStatSig(), 'Residuals:', peaks[i].GetResiduals(), '\n'
    print peaks[i].GetPeakStart().GetYErrLow(), peaks[i].GetPeakEnd().GetYErrLow()
print 'END OF PEAK SEARCH RESULTS'

c1.draw(h)
h1.setFill(1) # show background
c1.draw(h1)
peak,lreg = MakeGraph(peaks)

# always ignore the first point
for i in range(1,len(peak)):
  if (peaks[i].GetStatSig() <minSig): continue
  print 'Peak Number:', peaks[i].GetPeakNumber(), 'Peak Start:', peaks[i].GetPeakStart().x, 'Peak End:', peaks[i].GetPeakEnd().x, 'Peak Stat Sig:', peaks[i].GetStatSig(), 'Residuals:', peaks[i].GetResiduals(), '\n'
  print peaks[i].GetPeakStart().GetYErrLow(), peaks[i].GetPeakEnd().GetYErrLow()

  peakG=peak[i]
  peakG.setColor(Color.red)
  peakG.setSymbol(8)
  c1.draw(peakG)

  lreG=lreg[i]
  lreG.setColor(Color.red)
  lreG.setDrawLine(1)
  lreG.setSymbolSize(1)
  c1.draw(lreG)

  x = peaks[i].GetPeakEnd().GetX()
  y = peaks[i].GetPeakEnd().GetY()
  ssig = '#sigma = ' + str(round(peaks[i].GetStatSig(), 1))
  lab=HLabel(ssig,x+100,y+100, "USER")
  print ssig,x,y
  c1.add(lab)
  c1.update()

Run this code and you will see an image with identified peaks. It also calculate significance and the width. The output of this algorithm is shown in this image:

The identified peaks are shown in blue color.