# (Radiation) Pulse detection and height characterization when pulses are stacked

Edit: Title

Edit2: Clearer picture of low count rates.

Edit3: Added 2nd sentence for context.

I've been looking at some scope traces I took of radiation. This is a trace taken from a Geiger Muller (GM) tube. My DSP is really rusty and I can't even remember what to call some of this.

To start off, I've got a signal like the following, at low rates this is fairly simple to find the peak and identify the height. I just take the 503 - the peak.

Sometimes events stack up and at first I was content to just use the value of each point relative to the dotted line.

As the pulse rate gets higher this starts to become more of a problem. And my peak detection starts to miss more and more.

I can do a high pass filter to get rid of some of the noise and that helps a fair amount. However I'm having trouble finding the sweet spot between cleaning the signal, ringing at the peaks, and phase/amplitude shifting the signal. Below is a Butterworth lowpass.

I'd like to also be able to characterize the height based on the distance moved and I'm thinking I can find the high peaks and do some logic to try to make a reasonably accurate estimations.

My end goal is to get a histogram of the pulse heights and distances between pulses. The following is an example based on the peak point vs the purple dotted baseline shown in some of the other plots.

But I feel like I'm well into the weeds at this point and my question is this.

Does anyone have any suggestions on how to improve this or a different approach I can take? I have a vague recollection of some FIR (I think its FIR anyway) filters that can "transform" pulses like this but nothing is jumping out at me in my searches.

High level view. Note the "Bananas" level has so many micro pulses its making the large pulses consistently spaced.

Update

Deconvolution Attempt:

Following one of the suggestions below I tried to find pulses via deconvolution. I say attempt because I'm obviously missing something.

from scipy import signal

single_peak_l_edge = find_nearest_position(raw_dat["Time"], -0.000115)
single_peak_r_edge = find_nearest_position(raw_dat["Time"],  0.0002)
single_peak = raw_dat[single_peak_l_edge:single_peak_r_edge].copy().reset_index(drop=True)

decon, remainder = signal.deconvolve(raw_dat["Channel A"], single_peak["Channel A"])

fig, ax = plt.subplots(1, 1)
ax.plot(single_peak["Channel A"])
ax.set_title("Template")

fig, ax = plt.subplots(2, 1, sharex=True)

ax[0].plot(raw_dat["Channel A"])
ax[0].set_title("Origional")
#ax[0].set_xlim(0, 40e3)

ax[1].plot(decon)
ax[1].set_title("Deconvolved")


• Do the underlying pulses always have the same height and falloff rate? – endolith Mar 26 at 20:03
• "pileup mitigation in photomultiplier tube operation" (or avalanche photodiodes, or proportional counters) is an age-old problem. The things you'd like to read were written in the 1960's to 1990's. Pileup is a serious problem that you generally have to live with. From what I remember, at best you can detect and reject these events sometimes, but mostly you simply have to use lower count rates, either by reducing the events per second produced, or with granularity; using a larger number of detectors with smaller areas and lower count rate. See also my comment on (at)endolith's answer. – uhoh Mar 27 at 1:57
• Of course signal processing was done via analog circuits back then, and the necessary computing to analyze these pulses is readily available and cheap compared to then, but the fundamentals of pulse analysis and pulse shape variability due to charge collection statistics and other real effects will be just as much a problem for a digital solution as it is for an analog one. – uhoh Mar 27 at 1:59

Can you use deconvolution to convert these decaying impulses back into impulses?

Proof of concept:

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

impulses = np.zeros(150)

times =   [3,  20, 30, 40, 45, 50, 55, 80, 90]
heights = [-8, -7, -1, -9, -1, -2, -1, -8, -1]
impulses[times] = heights

b, a = signal.butter(1, 0.04)

filtered = signal.lfilter(b, a, impulses)
fig, ax = plt.subplots(3, 1, sharex=True)
ax[0].plot(filtered[:151])
ax[0].set_xlim(0, 100)
ax[0].set_title('Original')

impulse = -10*signal.unit_impulse(50)
template = signal.lfilter(b, a, impulse)
ax[1].plot(template)
ax[1].set_title('Template')

deconvolved, remainder = signal.deconvolve(filtered, template)
ax[2].plot(deconvolved)
ax[2].set_title('Deconvolved')


Then they are separated from each other and you can count their heights and events more easily.

Depends what the properties of the original bumps are, but they look like lowpass filtered impulses, which would mean they are all from convolution with the same impulse response, so this should work?

• This is what I was half remembering. I'm not sure if the "micro" pulses have a similar enough shape though. Am I correct in assuming you can deconvolve a second template as a separate operation to get separate set of impulses? I'll have to play around with this either way. Thanks! – user3219864 Mar 26 at 21:25
• @user3219864 are the micro pulses behind generated by the same process, just with smaller initial amplitude? – endolith Mar 26 at 22:56
• @user3219864 I'm guessing that because your voltage is positive you are collecting electrons produced in a vacuum (photomultiplier) a gas or liquid (proportional counter) or a solid (avalanche photodiode). In all these cases the pulse shapes can vary in a non-linear way due to space charge limitation and if it's a photomultiplier attached to some scintillator materials the ionization density variation due to particle type and energy can add another variable component to the pulse shape. – uhoh Mar 27 at 1:55
• @user3219864 Maybe you can add a little bit more to your question, indicating what the source of the signals are? If you are just counting individual photons of a constant wavelength with a PMT or avalanche counter the variation might be less than if you are using a radiation detector of some kind. See also my comment under your question. – uhoh Mar 27 at 1:56
• @uhoh Edited question also, but this is GM tube using a Cs137 source. I realize this is a bit like putting laser guidance on a hammer but its something I was playing around with and curious what might be possible with DSP. It sounds you you probably know this already but the typical method is to use a simple comparator to detect when the signal crosses a threshold. "Dead Time" is then mathematically added to account for time the detector is "dead" due to the recovery from an event. All measurement are based on the "count rate". Mostly this is me trying to practice some DSP. Maybe get insight. – user3219864 Mar 29 at 14:38

Your y-axis is in volts does the voltage increase overtime or does it fluctuate within a range? Without knowing more about the signals characteristics and what kind of frequency content you're dealing with it's hard to say..

You could possibly use a wavelet type approach to capture slope changes in your signal.

To get an idea of what I mean see this unanswered question I made that shows a plot of the continuous wavelet transform of a signal that has two "events" if you could call them that.

cwt plot example

• There is a DC level its always attempting to return to. I took traces at different field levels so each data set is "consistent". Yeah, I realize now I may have trimmed this screenshots a little too much trying to keep it concise. While working on this project I looked into this a bit. I'm pretty fascinated with them but they are a bit over my head at the moment (and I've probably exhausted most of the time I can spend on this little side project). If you'd like to play with the data I'd be happy to provide you with some. I'd be interested in what you came up with. – user3219864 Mar 26 at 21:39
• Added another pic to help with high-level view. – user3219864 Mar 26 at 22:07