[Tutor] Peaks detection

[Roeland Rengelink]

Marco wrote:

>Hey all, 
>i am tryin to write in python a peak detection algorithm to use on spectra taken from a data analyser in a lab.
>The spectrum comes in a text file with 2 columns, one with the channel number and the other one with number of counts for that channel.
>Making the calibration to the energy is not a problem, except for the fact that the algorithm _should_ recognize the peaks of the spectrum (those that have a fair gaussian form) and notate the centroid of the gaussian.
>
>The problem is: i can't find an algorithm except for a labview "procedure" on the net.
>Once found the "how" implementing it in python shouldn't be that hard :)
>  
>
>so the question is:
>Anyone has a source for an algorithm of that kind? Any suggestion?
>  
>
Hi Marco,

I'm only going to give you an outline of the solution, because this 
problem is a lot easier to solve in theory than in practice (code is 
untested).

Lets't define two functions:

def is_local_maximum(counts, current_channel):
    """Return True if counts in current_channel are larger than
       in neighbouring channels"""

    if (counts[current_channel] > counts[current_channel-1] and
        counts[current_channel] > counts[current_channel+1]):
        return True
    else:
        return False

def measure_peak(counts, channels, max_channel, width)
    """Return the mean (centroid) and stddev of a peak at max_channel"""

    mean = 0.0
    stdev = 0.0
    total_counts = 0.0
    for i in range(max_channel-width, max_channel+width):
        total_counts += counts[i]
        mean += counts[i]*channels[i]
    mean = mean/total_counts
    for i in range(max_channel-width, max_channel+width):
        stdev += counts[i]*(channels[i]-mean)**2
    stdev = math.sqrt(stdev/total_counts)
    return mean, stdev

To analyze the dataset we are going to assume that you know the 
background signal and the rms of the background signal (measuring them 
is left as an exercise)


def analyze(counts, channels, background, background_rms, min_sigma, 
window_width):
    minimum_signal = background+background_rms*min_sigma
    for i in range(window_width, len(channels)-window_width):
        if counts[i] < minimum_signal:
            # Ignore noise
            continue
        if is_local_maximum(counts, i):
           mean, stdev = measure_peak(counts, channels, i, window_width)
           print "Detected peak", counts[i], mean, stdev

Now for the caveats:

- I assume that there is going to be only one peak within the window 
around a given maximum
- I assume that the peaks are properly sampled (channel_width < 
width_of_peak/2)
- I assume that the window width >> width of the peak
- I assume that the total signal within the window is dominated by the peak
- I assume that the background is flat and the rms of the background is 
constant

These assumptions are generally not (all) true in practice

Hope this help,

Roeland