Geiger-Muller counter interfacing with Raspberry Pi.
| Article in progress |
When interfacing a Geiger-Muller device with a MCU or a Raspberry pi, there are mainly two options available :
- Use of a high level device, which transfers counter information digitally such as cps, cpm, tube type, stored memory or whatever information the device holds. Some devices may also be able to inform the host of a pulse detection using a low level signaling mode, and perform as a low level device. High level devices may be more or less sophisticated and allow for calibration, unit conversions, dosimetry, or changing settings depending on the tube model and their efficiency characteristics. High count compensation due to dead time may or may not be done, but any serious device should perform such processing. The compensation formulas are dependent on the designer knowledge of the matter and their proprietary choices.
- Use of a low level device. Such a device is an analog front end that supplies high voltage to the GM tube, and outputs detected events as pulses on a GPIO pin. The short pulse is usually stretched to a fixed time pulse, with sharp rising and falling edges. It should essentially output a square wave signal, that is easily processed by subsequent A/D stages that are not part of the device itself. Signal output is commonly done via GPIO header pins and/or 3.5 audio jacks. Most also feature a buzzer for direct audio pulse output.
This article is focused on the low level devices and their interfacing with a digital counting platform such as the Raspberry Pi, as it offers maximum processing flexibility.
Note that for high CPM processing, digital latency effects starts to be significative when using a non real-time OS platform due to operating system overhead. Part of this article will deal with mitigation of these unwanted effects.
GM Tube basics.
There are plenty of resources explaining the physics of the GM Tube. We won’t cover them in detail except for the important issue of dead time and how this skew the measured count rate at high levels of ionizing radiation (high count rates). This issue is critical as the device underestimates the count at these regimes. That is, the device is not linear.
The GM Tube is supplied a high enough voltage to work in the Geiger plateau, such as an ionization event triggers a Townsend avalanche and a significant charge is displaced between anode and cathode, and this is seen by the amplifier (either a transistor or an op-amp) as a voltage pulse. This charge displacement is equivalent to a lowering of the impedance of the GM Tube, which induces a voltage drop from the nominal high voltage, as the high voltage source is far from ideal, it can’t source much current and has added resistors. The GM tube and resistors on the high voltage path perform as a voltage divider (with a parallel path to ground, the GM tube – with its impedance lowering during the Townsend avalanche process)
This voltage drop is conditioned to a low voltage (with reference to ground), buffered, and the pulse is stretched to several hundreds of microseconds to a couple of milliseconds, depending on the design and/or tuning.
Pulse stretching usually involves stages such as Schmitt triggers and 555 timers
Dead time from the analog front end.
Usually the analog front end circuitry processing exhibits a non-paralyzable dead-time characteristic.
dead time : it means that after a pulse is detected by this front end, any subsequent pulse is not detected during the interval spanning from the registered pulse to the expiration of the dead-time (the pulse stretched “HIGH’ state time, set usually by an analog timer such as a 555.)
non-paralyzable : it means that any townsend avalanche event occuring during that dead-time interval. does not re-trigger the dead time, or at least, it should not. If it does, it would be a poor design choice as dead-time from any processing stage should be minimized as much as possible, and a paralyzable dead-time may extend this duration.
GM Tube intrinsic dead time.
As for the GM tube Townsend avalanche process, it has an intrinsic dead time characteristic that is short but still significative. After a sucessful ionizing event detection, one that triggers a full Townsend avalanche, the tube is unable to detect any other event during the dead time. This dead time is usually in the order of 80 to 200 µs, and is mostly paralyzable : Any event ocurring during that dead time, resets that intrinsic dead time to some extent. For sake of simplicity, we will treat it as a fully paralyzable dead time. In reality, most of the literature agrees that this intrinsic dead time falls in between a paralyzable and non paralyzable model.
We make the supposition that the dead time is a function of the time separation of two subsequent (or more) ionization events, with the limit reaching the nominal dead time as event separation in the time domain goes toward infinity.
The precise model of a Townsend avalanche and its practical implementation in a commercial GM tube is outside of the scope of this article.
Thus it follows that the dead-time nominal values from the GM Tube and those of the analog pulse stretch front end are important parameters for count rate compensation.
Analog frontend deadtime measurement is quite straightforward in a low count (low background radiation) environment, and can be done by plugging the GM meter jack output to a sound card set at the highest sampling capability, ideally 96ks/s or more. Or better, by using an oscilloscope.
In our case, the dead time was measured at around 2.2ms average on this commercial device.
As for the intrinsic dead time of the GM Tube, some are reported in the datasheets. Determination of that dead time is not an easy task, in part because the ionization events exhibit a probabilistic distribution (modeled as a Poisson process) – That is, it is hard to control the precise moment a ionization event will trigger a Townsend avalanche.
Furthermore, GM Tubes have an efficiency parameter: Of all the gamma flux that intercepts the section area of the tube (usually the tube length times the tube diameter, and taking into account the mass absorption characteristics of the tube wall due to limb effects if one wants a really precise characterization), not all gamma flux will interact with the tube gas electrons, the main factors to take into account are the finite – and small – tube volume, low gas density, and type of gas used as it influences ionization behaviour. Note that GM tubes are sensitive to x-ray (with quite low efficiency) gamma rays, and charged particles such as alpha, beta+ and beta-. Non charged particles such as neutrons are detected indirectly, by a nuclear process that give rise to charged particles when a neutron is captured or absorbed by the GM tube gas. Neutron sensitive GM tubes may use BF3 gas (Boron Tri-Fluoride)
GM Tube efficiency
Common beta/gamma small glass tubes exhibit an efficiency of around 0.03
Moreover, this efficiency varies with the gamma ray energy. It is well known that GM Tubes under-perform in the X-ray spectrum. Flattening of the energy/efficiency response curve can be done by using filters (such as thin metal sheets) that partially occult the tube, mainly for gamma counting.
As for tubes that exhibit a better sensitivity in the alpha and beta detection mode, they usually are made with walls that have a lower extinction coefficient for alpha or beta particle, such as mica.
Remember also that GM tubes are not proportional counters, scintillators, or particles chambers. They do not convey information about the source of the ionization event energy, as the Townsend avalanche “saturates” the detector.
True count rate estimation based on GM tube dead time and analog front end dead time
Based on the research of Muller in “Dead Time problems” – Nuclear instruments and methods 112th issue (1973), (1) The model we will use is stated on page 56 (d) equation (56), as it accounts for two dead times in series, one paralyzable from the GM tube, and one non-paralyzable from the analog front end (pulse stretcher)
we can see that the analog dead time due to pulse shaping into a square wave is the limiting factor as for high counting rates. Thus, it is preferable to tune, in the design phase, the 555 timing to achieve the shortest pulse length that does not trigger spurious counts by the A/D stage. A/D stages typically register the rising edge of the pulse. As the amplitude of the pulse is not of any particular interest, a high frequency, low bit resolution A/D design is preferable, such as a sigma/delta ADC, (ADC considerations go beyond our scope, as we are limited by the BCM GPIO hardware on the Raspberry Pi) and good noise immunity practices, such as using shielding for the whole assembly and a shielded signal cable grounded at one end, or twisted pair, while maintaining the shortest possible signal path from the GM front end to the ADC. This is to ensure good signal integrity and low stray loop inductance. A very short pulse could induce ringing and be registered as multiple events when the level lingers between the low/high logic boundary. However, there are diminishing returns for high count precision when the analog front end dead time approaches the intrinsic GM tube dead time, that is in the order or 80 to 200 µs. These time constants are well within the 555 timer capabilities, and proper design should take care of ringing artifacts. As for the hardware limitations of Raspberry pi pulse sensing, they are in the tens of ns range. Which means that lowering analog front-end dead time for ADC sensing to around 100 µs should not be a hardware limiting factor. As for the software GPIO / ISR library, it is recommended to use one that minimizes delay. Some libraries are implemented as linux daemons that may allow better management of high count rates. Best is to compare performance with a signal generator between the standard Python Raspberry Pi GPIO libraries that perform better with high pulse frequency (in the order of 10 kHz).
Python code processing could be the bottleneck, if care is not taken to benchmark deque() performance at 10 kHz pulse rate, particularly if the deque() stores the pulse timing with calls to Python time libraries.
Some GPIO pulse processing libraries store pulse timestamp information, in this case this overhead could be discarded.
Also, registering the falling edge may be useful to assess proper operation of the GM analog front end. Absence of a falling edge in a timely manner before the next rising edge could signify that the signal is stuck in a high state, and should display a malfunction and/or high count warning.
It is also preferable to use a micro controller or full fledged miniature computer such as a Raspberry pi with adequate processing power, which translates into a fast CPU clock and more than one core (for computers) to decrease ISR (interrupt service request) burden to a minimum in high count environment.
In the case of a Raspberry Pi, and when using Python, the following guidelines should be followed :
- Minimum amount of code in the ISR.
- Ideally it should use a deque() for pulse registering, simply appending the pulse to the deque in the ISR, and exiting the ISR.
- Assessing the need for specialized GPIO libraries for pulse counting, such as those that involve a daemon, when pulse timestamping with the highest precision is a design requirement. In that case however, inter-process communication (IPC) delay is a factor to take into account for overall responsiveness.
- A separate thread on another core should perform post processing such as logging into the file-system or a database, or count rate compensation calculations.
- It is preferable to use a low level language such as C++ instead of Python, and precise benchmarking should be done using a function generator generating pulses with a comparable duty time to the GM pulse shaper backend, with increasing CPM frequency, to assess the digital latency induced and ISR responsiveness in high CPM environments. This way, the influence of the A/D backend and code performance can be precisely factored in for final count up-rating & device calibration.
Pulse time-stamping
In our design the deque() stores a timestamp for each pulse. Pulse time-stamping is a niche requirement and may be of use to triangulate the source of a sudden burst of gamma energy, although air attenuation coefficient factor coupled with poor efficiency of GM tubes, would render such an endeavour tricky. As for alpha and beta particle radiation, their detection on fixed position sensors for environmental monitoring in weather stations is dependent on wind and atmospheric currents drift that carry contaminated dust and fallout, which operate on timescales large enough to not require precise time stamping.
For EAS (extensive air shower) research arising from high energy cosmic particles, as well as study of dark lightning generated terrestiral gamma-ray flashes (TGF) a scintillation detector would be a better choice, as they have a better quantum yield.
A GPS module is a good investment in a project of this kind as it allows not only Geo tagging of events, but also precise timestamping due to inherent time synchronisation features of GPS.
Alternatively, low quality synchronisation of nodes that perform event timestamping may use NTP or higher precision NTP protocols. In any case, precise node time synchronisation using NTP requires symmetric network packet processing (no asymmetric routing – this creates different propagation delays upstream and downstream). These propagation delay uncertainties increase substantially when the time source is several router hops away, and also depend on the network traffic load induced delay.
If time synchronisation is performed through air (such as using Lora) all radio induced propagation delays and bitrate induced delays have to be factored in.
Count up-rating using two dead times in series model with time constants t1 and t2.
In our project, we will use Muller’s derivation (1) p56. (d)
$$ R = \frac{\rho}{(1-\alpha)x} + e^{\alpha x} $$
$$ x = \rho t_{2} $$
$$ \alpha = \frac{t_{1}}{t_{2}} $$
$$ \rho $$
being the true count rate estimation, and
$$ R $$
being the measured count
$$ t_{1} $$
being the (paralyzable) dead time of the GM Tube, and
$$ t_{2} $$
being the (non-paralyzable) dead time of the analog frontend pulse shaper.
Let’s introduce the other well known models accounting for a single dead time system:
the non paralyzable model :
$$ R = \frac{\rho}{1 + t\rho} $$
the paralyzable model :
$$ R = \rho e^{(-\rho t)} $$
Since the unknown is the corrected count (rho), we need to use the inverse function of these models, regardless of the model, compound, paralyzable or non paralyzable.
The paralyzable function inverse expression requires the use of the W0 and W1 Lambert function, Math helpers in Python such as scipy allow straightforward calculation of the Lambert W0 and W1 branches, albeit with some computational burden.
The compound t1 and t2 in series requires numerical methods such as the secant method. Which would only increase the computational burden.
In the case of a Raspberry Pi, since RAM and storage are not an issue, and the problem is not multivariate since t1 and t2 are constants. We advise to compute the functions models, and use a reverse table lookup for fast determination of the corrected count. Scipy propose linear and higher order interpolation mechanisms, which would have a lower computational burden than root finding.
Calibration and experimental confirmation of the count correction model.
Disclaimer : For this part, access to a medium activity point source > 100 kBq, with low uncertainty is preferable. A low activity source would not push the GM counter at count rates where there is significant deviation from the true rate, and thus not give enough data to properly test the models. Depending on your location, point sources are exempt of declaration at different ceilings. To complicate things further, several regulations may overlap such as national / European Union, and providers of sources may or may not be available in your country. a source in the 100 kBq pretty much requires a license anywhere in the world. One option may be to use the source of a third party at the site of the third party licensed metrology / research lab. As always, exposure to a medium activity source may be harmful depending on exposure time and proper handling and shielding precautions are required.
Besides the point source, a lead plated rectangular cross section channel spanning from the point source to the GM Tube is preferable, such a device function is not to serve as a collimator, but rather to prevent reflections of gamma rays leading to constructive / destructive interference, as some papers suggest. It seems quite remarkable however that such effects interfere significantly with the measurement, given the very low refractive index of most materials in the gamma spectrum. They would however shield the detector somewhat from background radiation and other sources that may be found in a lab at relative proximity.
Note that a Cs137 nucleide decays either to a stable 56Ba137 nucleide through Beta- decay, with a probability of 5.4%, or to an excited (56Ba137m) state, also through Beta- decay, with a probability of 94.6%. When reverting to the ground state 56Ba137, some of the 56Ba137m nucleides emit a 0.6617 MeV gamma photon. Of all Cs137 decays that form the measure of the total activity of the sample supplied by the manufacturer, 85.1 % yield gamma photons. This has to be factored in, besides manufacturer supplied source percent uncertainity, and exponential law derating of the activity of the sample, using the supplied date of manufacture, and the Cs137 half-life of 30.05 years.
Beta- radiation has a high linear attenuation coefficient in air, and would skew the measurements when the jig is very close to the tube through their contribution in the final measure. A standard food grade aluminium foil is sufficient to filter them out to a negligible level, while leaving almost all gamma ray photons unscathed. Nevertheless we have factored in attenuation from air and aluminium through air and aluminium for both Beta- and gamma based on litterature and NIST database data.
The source is then placed on a jig powered by a linear actuator, so that it moves freely inside the channel along the x axis. the center of the point source disc should be on the x axis and intercept the center of mass of the GM tube. The x axis should be perpendicular to the point source disc and GM tube.
The raspberry Pi platforms registers pulse and operates the linear actuator. A model with position feedback is required for accurate jig position determination, after manual position calibration, or alternatively, a time of flight (TOF) / LIDAR sensor module should be used. This solution should also be subjected to calibration beforehand.
The GM calibration helper for a Cs-137 point source Python script is in alpha stage of development but the bulk of the work is done.
The workflow of the script starts with the following parameter inputs :
- Point source activity, supplied in either Bq or Ci units.
- The activity_unit enumeration is used to set up the unit of activity
- The radioisotope used is hard-coded to Cs137, and it’s decay scheme is factored in to take into account only Gamma activity.
- Decay % that give rise to gamma photons
- Cs137 half-life as a constant
- Source date of manufacture
- Aluminium mass extinction coefficients for Beta and gamma radiation (using Cs137 decay energies)
- Air mass extinction coefficients for Beta and gamma radiation (using Cs137 decay energies)
- Distance range (min,max) of the source relative to the GM tube center of mass on the jig y-axis
- GM Tube geometry : length and diameter.
- GM tube intrinsic dead time estimation (or from datasheet)
- GM tube efficiency estimation (or from datasheet)
On that basis, we will first calculate the mean path length from the point source to the tube as a function of y axis distance. For this we will use basic trigonometry and integrate over the tube length, all the ray paths from the source. This mean path length will be used to factor in the linear attenuation coefficient of air for gamma radiation at 0.611 MeV energy. Beta radiation is assumed to be filtered up to an insignificant amount, but we can still calculate the attenuation based on aluminium mass extinction coefficient and foil thickness to check for filtering efficiency of beta radiation.
The compound attenuation of air plus aluminium filter is calculated for the gamma radiation.
The radiant flux from the point source received by the GM tube is then calculated. The standard law
$$ Flux =\frac{P}{4 \pi r^{2}} $$
assumes a spherical irradiated surface. since the cross section of the GM Tube is a plane, and taking into account that $$ GMtubewidth \ll GMtubelength $$, we will perform a single integration along the tube cylinder axis to get the corrected number of photons $$ R_{net} $$ crossing the tube section, instead of a double integration.
$$ R_{totvac} \simeq 2\times\int_{0}^{\frac{GMtubelength}{2}}\frac{P_{net}}{4 \pi (x^2 + r^2)} cos(arctan(\frac{x}{r}))dx \hspace{2mm} \times GMtubewidth $$
Which simplifies to :
$$ R_{totvac} \simeq 2\times\frac{P_{net}\frac{GMtubelength}{2}}{4\pi r \sqrt{(\frac{GMtubelength}{2})^{2} + r^{2}}} \times GMtubewidth $$
Also,
$$ P_{net} $$ is the gamma source activity. It has to take into account the percent of decays giving rise to gamma photons, and the reduction of activity since the source date of manufacture.
$$ P_{net} = 0.851 P_{nom} e^{-\frac{ln(2)}{t_{1/2}} (t_{1} – t_{0})} $$
$$ t_{1} $$ is the time of the calibration
$$ t_{0} $$ is the time of source manufacture.
and, accounting for attenuation :
$$ R_{tot} \simeq 2\times att_{\gamma}\frac{P_{net}\frac{GMtubelength}{2}}{4\pi r \sqrt{(\frac{GMtubelength}{2})^{2} + r^{2}}} \times GMtubewidth $$
To calculate the air attenuation of the gamma photons, we will use an approximation that has an expression in terms of elementary functions.
A gamma photon path length ‘l’ from the point source to any point on the x axis (axis of the GM tube represented as a cylinder) can be approximated as the length of the hypotenuse (neglecting paths that strike the cylinder off axis) :
$$ l = \sqrt{x^{2} + r^{2}} $$
The average path length is :
$$ l_{avg} = {\frac{2}{GM tube length}} \int_{0}^{\frac{GM tube length}{2}} \sqrt{x^2 + r^2}dx $$
Which yields the following expression :
$$ l_{avg} = \frac{{r^{2} \operatorname{arsinh}\left(\frac{GMtubelength}{2r}\right)}}{GMtubelength} + \frac{\sqrt{({\frac{GMtubelength}{2}})^{2} + r^{2}}}{2} $$
The attenuation factor applied to the flux :
$$ att_{\gamma}\simeq e^{-\mu_{air} l_{avg} -\mu_{Al}w} $$
$$ w $$ being the Aluminium foil thickness
$$ \mu_{air} $$ being the linear attenuation coefficient of air at 0.661 MeV gamma energy.
$$ \mu_{Al} $$ being the linear attenuation coefficient of Aluminium at 0.661 MeV gamma energy.
A small fraction of these photons will trigger an ionization of the GM_tube gas medium. These will form the measured tube event count per second. $$ R_{net} $$ At low count rates, dead time effects are negligible, and provided sufficiently long measurement times, the tube efficiency alpha can be determined :
$$ \alpha = \frac{R_{net}} {R_{tot}} $$
Resources
Note : For up to date code, please check https://www.github.com/rodv92/GMobile
This is the calibration helper Python script draft.
import scipy as sp
import math as m
import numpy as np
import time as time
from datetime import datetime
import matplotlib.pyplot as plt
from enum import Enum
# using S.I units (unless specified otherwise in comment immediately following
# declaration)
class activity_unit(Enum):
BQ = 1 #Bq
CI = 2 #Ci
activity_unit = Enum('activity_unit', ['BQ','CI'])
# assuming a Cs-137 source, and a collimator made of thick lead of the same width as the GM tube and the same height as the tube, with the point source at the e# ntrance of the channel
P_nom = 370e-3 # source activity (unit determined by used_activity_unit)
used_activity_unit = activity_unit.BQ
# it's better to use a high activity source in order to drive the count at very high levels for dead-time effects to be come noticeable
# the lower the analog frontend dead-time that results in reliable A/D conversion, the higher the activity of the source is required.
gamma_yield = 0.851 # ratio of decays that give rise to gamma photons
source_date_of_manufacture = "2022-12-01T19:00:00Z"
dt_source_dom = datetime.fromisoformat(source_date_of_manufacture)
dt_now = datetime.now()
source_age = (dt_now - dt_source_dom).total_seconds()
half_life = 30.05 # Cs137 half-life in years
half_life = half_life*365*24*60*60 # Cs137 half-life converted to seconds
µ1 = 0.000103 # linear attenuation coefficient of air for the gamma energy of 667 kEv - emmited by Ba137m
#µ2 = 0.1 # linear attenuation coefficient of beta shield for the gamma energy of 667 kEv. We assume that the beta shield filters beta particles to a negligible amount.
Al_filter_thickness = 0.001 # cm 0.001 cm = 10µm standard food grade aluminium foil thickness
max_distance = 0.3
distance = 0.15 # distance from point source to center of detector. (meters)
min_distance = 0.05
# GM tube axis is normal to line from point source to GM tube center.
GM_tube_length = 0.075 # glass envelope length only (meters)
GM_tube_diameter = 0.01 # glass envelope diameter (meters)
GM_tube_alpha = 0.0319 # tube efficiency - not all incoming photons trigger an ionization event.
alpha = GM_tube_alpha
GM_tube_detection_cross_section_area = GM_tube_length*GM_tube_diameter
GMT_dcsa = GM_tube_detection_cross_section_area # variable alias
GM_tube_dead_time = 80e-6 # tube intrinsic (charge drift time) dead time
GMT_det = GM_tube_dead_time
if (used_activity_unit == activity_unit.BQ):
pass
elif (used_activity_unit == activity_unit.CI):
P_nom *= 3.7e10 # convert to Bq
#derate activity based on % decay gamma yield and source age.
P_net = gamma_yield*P_nom*exp(-(m.log(2)/half_life)*source_age)
#calculate mean path length of gamma photons reaching the tube.
mean_path = (1/2)*(distance**2 + (GM_tube_length/2)**2)**(1/2) + (distance**2)*m.arcsinh*(GM_tube_length/(2*distance))/GM_tube_length
#calculate attenuation from linear attenuation coefficient
gamma_att_air = m.exp(-µ1*mean_path)
Al_density = 2.7 # g/cm3
#Al_mass_att_beta_cs137 = 15.1 #cm2/mg
Al_mass_att_beta_cs137 = 15.1e3 #cm2/g beta+/- attenuation for Cs137 emitter https://doi.org/10.1016/j.anucene.2013.07.023
# not used in subsequent calibration formulas, assuming that beta is filtered to an insignificant amount. standard food grade aluminium foil of 10µm thickness gives an attenuation factor in the order of 1e-18
µ2 = Al_mass_att_beta_cs137*Al_density #cm-1
beta_att_Al = m.exp(-Al_filter_thickness*µ2)
# To check the effectiveness of beta filtering.
# not used in subsequent calibration formulas, assuming that beta is filtered to an insignificant amount. standard food grade aluminium foil of 10µm thickness gives an attenuation factor in the order of 1e-18
Al_mass_att_gamma_Ba137m = 7.484e-2 #cm2/g gamma attenuation for Aluminium at 667 kEv
µ3 = Al_mass_att_gamma_Ba137m*Al_density
gamma_att_Al = m.exp(-Al_filter_thickness*µ3)
#calculate total gamma attenuation from air and Al filter.
gamma_att_total = gamma_att_air*gamma_att_Al
import GMobile as GM
# CALIBRATION step 0
# Estimates the analog front end dead-time by timing the pulse width while the GM Tube is exposed to background radiation
background_measure_pulsewidth_total_pulses = 60 # time to spend in seconds measuring pulse width
background_measure_pulsewidth_max_fails = 20 # time to spend in seconds measuring pulse width
rise_timeout_ms = 20000
fall_timeout_ms = 20
(pulsewidth,stdev) = GM.measurePulseWidth(background_measure_pulsewidth_total_pulses,background_measure_pulsewidth_max_fails,rise_timeout_ms,fall_timeout_ms)
# CALIBRATION step 1
#Measure the background CPM with the collimating assembly, but the source removed and far away.
# call main.py and average CPM over specified background_acquire_time in seconds
background_acquire_time = 600 # time to spend in seconds acquiring background radiation levels after first 60 sec of acquisition.
chars = None
chars = input("Step 1 - acquiring background radiation cpm during " + background_acquire_time + " seconds. Please put the source as far away as possible. press ENTER to start")
while chars is None:
chars = input("Step 1 - acquiring background radiation cpm during " + background_acquire_time + " seconds. Please put the source as far away as possible. press ENTER to start")
GM.SetupGPIOEventDetect() # sets up the GPIO event callback
s = 0
cpm_sum = 0
while (s < background_acquire_time):
cpm = GM.process_events(False,False) # This call should take exactly one second.
if (cpm != -1):
cpm_sum += cpm
s += 1
cpm_background = cpm_sum/background_acquire_time
def model1_estimated_GM_CPM(true_count,t1,t2): # Muller's serial t1 paralyzable dead time followed by t2 non paralyzable dead time model
alpha = t1/t2
x = true_count*t2
corrected_cpm_1 = true_count/((1-alpha)*x + m.exp(alpha*x))
return corrected_cpm_1
def model2_estimated_GM_CPM(true_count,t1):
corrected_cpm_2 = true_count*m.exp(-true_count*t1)
return corrected_cpm_2
def model3_estimated_GM_CPM(true_count,t2):
corrected_cpm_3 = true_count/(1+true_count*t2)
return corrected_cpm_3
def movejig(position):
#TODO : linear actuator positioning code
err = 0
return err # err = 0 : actuation OK
def efficiency_step(distance=0.15,movestep=0.005,efficiency_placing_time=600,efficiency_stab_time=60,last_secs_stab_time=60,min_cpm_efficiency_cal=8*cpm_background,max_cpm_efficiency_cal=16*cpm_background):
# CALIBRATION step 2
print("Step 2 : automated GM tube efficiency calculation")
s = 0
s_stab = 0
cpm_stab = []
while (s < efficiency_placing_time):
#print("Step 2.1 - efficiency computation: you have " + efficiency_placing_time + "seconds to put the source at a distance to obtain a reading between " + min_cpm_efficiency_cal + " and " + max_cpm_efficiency_cal + " cpm. press ENTER when in range")
#print("the timer will start counting down after the first 60 seconds have elapsed - it will then exit the step if cpm is stabilized in range for a whole " + efficiency_stab_time + " secs, and get the cpm average for the last " + last_secs_stab_time + " secs. Do not move the source during that time")
cpm = GM.process_events(False,False) # This call should take exactly one second.
if (cpm != -1): # first 60 seconds have elapsed.
cpm_sum += cpm
s += 1
if(cpm > min_cpm_efficiency_cal and cpm < max_cpm_efficiency_cal): # in range.
s_stab += 1
cpm_stab.append(cpm)
elif(cpm > max_cpm_efficiency_cal): # out of range high. reset stabilized countrate time counter
distance += movestep
movejig(distance)
s_stab = 0
cpm_stab = []
elif(cpm < min_cpm_efficiency_cal): # out of range high. reset stabilized countrate time counter
distance -= movestep
movejig(distance)
s_stab = 0
cpm_stab = []
print("cpm:\t" + cpm)
print("stabilized_time:\t" + s_stab)
print("distance:\t" + distance)
if s_stab >= efficiency_stab_time:
#distance = input("Please input distance in meters from GM_tube to source at stabilized reading")
cpm_stab = cpm_stab[-last_secs_stab_time:]
cpm_stab_avg = sum(cpm_stab)/last_secs_stab_time
break
cpm_efficiency_calc = cpm_stab_avg - cpm_background
#calculate mean path length of gamma photons reaching the tube.
mean_path = (1/2)*(distance**2 + (GM_tube_length/2)**2)**(1/2) + (distance**2)*m.arcsinh*(GM_tube_length/(2*distance))/GM_tube_length
#calculate attenuation from linear attenuation coefficient
gamma_att_air = m.exp(-µ1*mean_path)
#calculate total gamma attenuation from air and Al filter.
gamma_att_total = gamma_att_air*gamma_att_Al
flux = gamma_att_total*2*GM_tube_width*(P_net*(GM_tube_length/2))/(4*m.pi*distance*((GM_tubelength/2)**2 + distance**2)**(1/2))
# gamma flux in photons.s^-1 crossing the GM tube. Accounting for planar cross section of GM Tube (instead of solid angle) and attenuation from air and beta filter of gamma photons.
theoretical_cpm = 60*flux # assuming GM tube efficiency of 1, All photons would give rise to ionization events inside the GM tube
if(s != efficiency_placing_time):
efficiency = cpm_efficiency_calc/theoretical_cpm
print("computed efficiency:\t" + efficiency)
return((efficiency,distance))
else:
print("efficiency calibration failed.")
return((-1,distance))
# Compute efficiency of detection at a count rate sufficiently high enough above background but not as high as dead time effects become significant.
efficiency_placing_time = 600
efficiency_stab_time = 180
last_secs_stab_time = 60
min_cpm_efficiency_cal = 8*cpm_background
max_cpm_efficiency_cal = 16*cpm_background
chars = None
chars = input("Step 2.0 - efficiency computation: please put the source at " + distance + "cm from the GM tube, normal to the tube, withint the collimator. press ENTER to start")
movejig(distance) # initial position
(efficiency,distance) = efficiency_step(distance) # gets efficiency, -1 if failed, and jig to source distance
while (efficiency == -1):
print("efficiency calculation step failed. repeating")
(efficiency,distance) = efficiency_step(distance)
#STEP 3 : record countrate while stepping the source jig towards the GM tube
step3_wait_time = 120
step3_last_seconds_measure = 60
step3_last_seconds_measure = min(120,step3_last_seconds_measure) # ensure the total seconds we sample is lower than step3_wait_time
x_distance = np.arange(max_distance,min_distance,-0.005)
y_cpm_m = np.empty(len(x_distance)) # array of average of count rate for each measurement sampling
y_cpm_m[:] = np.nan
y_std_m = np.empty(len(x_distance)) # array of standard deviation of count rate for each measurement sampling
y_std_m[:] = np.nan
y_cpm_t = np.empty(len(x_distance)) # array of theoretical cpms derated with GM tube efficiency
y_cpm_t[:] = np.nan
y_cpm_tm = np.empty(len(x_distance)) # array of theoretical cpms derated with GM tube efficiency and dead time effects
y_cpm_tm[:] = np.nan
distance_cpm_avg_m = np.stack([x_distance,y_cpm_m],axis=1) # tabular data for cpm (average) as a function of source/GM tube distance
distance_cpm_std_m = np.stack([x_distance,y_std_m],axis=1) # tabular data for cpm (std dev) as a function of source/GM tube distance
distance_cpm_t = np.stack([x_distance,y_cpm_t],axis=1) # tabular data for cpm (theoretical), derated by GM tube efficiency as a function of distance
distance_cpm_tm = np.stack([x_distance,y_cpm_tm],axis=1) # tabular data for cpm (theoretical), derated by GM tube efficiency and accounting for dead time effects
distance = max_distance # retract actuator to minimum to get longest source to GM tube distance.
if not (movejig(distance)): # no actuation error
idx = 0
while(distance > min_distance):
if(movejig(distance)):
break # actuation error, break loop
distance -= 0.05
cpm_stab = []
#TODO : reset deque() in GMobile after jig move, to get rid of the inertia induced by the sliding window sampling
while(s < step3_wait_time):
if(s >= (step3_wait_time - step3_last_seconds_measure)):
cpm_stab.append(GM.process_events(False,False)) # This call should take exactly one second.
else:
time.sleep(1)
cpm_avg = np.average(cpm_stab)
cpm_std = np.std(cpm_stab)
distance_cpm_avg_m[idx][1] = cpm_avg
distance_cpm_std_m[idx][1] = cpm_std
#calculate mean path length of gamma photons reaching the tube.
mean_path = (1/2)*(distance**2 + (GM_tube_length/2)**2)**(1/2) + (distance**2)*m.arcsinh*(GM_tube_length/(2*distance))/GM_tube_length
#calculate attenuation from linear attenuation coefficient
gamma_att_air = m.exp(-µ1*mean_path)
#calculate total gamma attenuation from air and Al filter.
gamma_att_total = gamma_att_air*gamma_att_Al
flux = gamma_att_total*2*GM_tube_width*(P_net*(GM_tube_length/2))/(4*m.pi*distance*((GM_tubelength/2)**2 + distance**2)**(1/2))
# gamma flux in photons.s^-1 crossing the GM tube. Accounting for planar cross section of GM Tube (instead of solid angle) and attenuation from air and beta filter of gamma photons.
theoretical_cpm_eff = 60*flux*efficiency # theoretical cpm (derated with GM tube efficiency estimated in step 2.0)
distance_cpm_t[idx][1] = theoretical_cpm_eff
distance_cpm_tm[idx][1] = model1_estimated_GM_CPM(theoretical_cpm_eff) # theoretical cpm from above with dead time compensation
def compare_cpm_measured_theoretical(theoretical,measured):
devsum = 0
devratiosum = 0
if (len(theoretical) != len(measured)):
return (-1,-1)
for idx in range(0,len(theoretical)):
devratiosum += abs((measured[idx][1] - theoretical[idx][1])/measured[idx][1])
devsum = (measured[idx][1] - theoretical[idx][1])**2
mape = devratiosum/len(theoretical)
return (devsum,mape)
def plotcurves(curve_x, curve_y1, curve_y2, color_curve_1, color_curve_2, label_curve_1, label_curve_2):
plt.plot(curve_x, curve_y1, color_curve_1, label=label_curve_1)
plt.plot(curve_x, curve_y2, color_curve_2, label=label_curve_2)
plt.xlabel('distance (m)')
plt.ylabel('count rate (cpm)')
plt.legend()
plt.grid()
plt.show()
print(compare_cpm_measured_theoretical(distance_cpm_tm,distance_cpm_avg_m))
plotcurves(distance_cpm_t[:,0],distance_cpm_avg_m[:,1],distance_cpm_tm[:,1],'r-', 'g-', "measured","theoretical_compensated")
This is the main GM counter interfacing script, with Mariadb logging, individual pulse timestamping Geotagging, and Lora telemetry (work in progress)
# coding: utf-8
import RPi.GPIO as GPIO
import scipy as sp
import signal
import sys
import time
import datetime
from collections import deque
import struct
import statistics
# Module Imports
import mariadb
#from mariadb.connector.aio import connect
import sys
import re
# Lora Imports
from SX127x.LoRa import *
#from SX127x.LoRaArgumentParser import LoRaArgumentParser
from SX127x.board_config import BOARD
def signal_handler(sig, frame):
GPIO.cleanup()
sys.exit(0)
def pulse_detected_callback(channel):
global pulse_events
pulse_events.append(time.time() - start_time_epoch)
def dms2dec(dms_str):
dms_str = re.sub(r'\s', '', dms_str)
sign = -1 if re.search('[swSW]', dms_str) else 1
numbers = [*filter(len, re.split('\D+', dms_str, maxsplit=4))]
degree = numbers[0]
minute = numbers[1] if len(numbers) >= 2 else '0'
second = numbers[2] if len(numbers) >= 3 else '0'
frac_seconds = numbers[3] if len(numbers) >= 4 else '0'
second += "." + frac_seconds
return sign * (int(degree) + float(minute) / 60 + float(second) / 3600)
# Connect to MariaDB Platform
try:
conn = mariadb.connect(
#conn = connect(
user="gmobile_user",
password="gmobile_passwd",
host="127.0.0.1",
port=3306,
database="gmobile"
)
except mariadb.Error as e:
print(f"Error connecting to MariaDB Platform: {e}")
sys.exit(1)
# Get Cursor
cur = conn.cursor()
start_time_epoch = time.time()
pulse_events = deque()
cps = 0
cpm = 0
#process_pulses = True
#paralyzable model
# m = n*exp(-nt)
# using lambert(w)
#y = w * exp(w)
#-nt = w
#y = -nt * exp(-nt)
#-nt = lambert(y)
#n = -lambert(y)/t
#scipy.special.lambertw(z, k=0, tol=1e-8)
#while True:
# GPIO.wait_for_edge(PULSE_GPIO, GPIO.FALLING)
# print("Button pressed!")
# sp.special.lambertw(z, k=0, tol=1e-8)
PULSE_GPIO = 16
lat = dms2dec("51°24'04.30\"N")
long = dms2dec("30°02'50.70\"E")
def setupGPIOEventDetect():
GPIO.setmode(GPIO.BOARD)
GPIO.setup(PULSE_GPIO, GPIO.IN)
GPIO.add_event_detect(PULSE_GPIO, GPIO.RISING, callback=pulse_detected_callback)
def removeGPIOEventDetect():
GPIO.remove_event_detect(PULSE_GPIO)
def measurePulseWidth(total_pulses,max_fails=20,rise_timeout_ms=20000,fall_timeout_ms=20):
removeGPIOEventDetect()
s = 0
fail = 0
pulsewidths = []
while(s < total_pulses and fail < max_fails):
channel = GPIO.wait_for_edge(PULSE_GPIO,GPIO.RISING,timeout=rise_timeout_ms) # assuming there is at least one pulse registered every rise_timeout_ms
if (channel is None):
fail += 1
continue
pulsestart = time.time_ns()
channel = GPIO.wait_for_edge(PULSE_GPIO,GPIO.FALLING,timeout=fall_timeout_ms) # assuming the pulse width is less fall_timeout_ms
if (channel is None):
fail += 1
continue
pulsewidths.append((time.time_ns() - pulsestart)/1000) # save pulsewidth value in µs
s+= 1 # sucessful pulsewidth measure
#this pulse timing method may be problematic if the real time separation between RISING and FALLING edge is shorter than execution time between the two
#wait_for_edge calls. thread should be set to high priority. Overall, dead time will be overestimated.
setupGPIOEventDetect()
return (pulsewidths[int(len(pulsewidths)/2)],statistics.pstdev(pulsewidths)) # returns mean pulsewidth value in µs and standard deviation.
# LORA INIT
#BOARD.setup()
#BOARD.reset()
##parser = LoRaArgumentParser("Lora tester")
"""
class mylora(LoRa):
global current_data
def __init__(self, verbose=False):
super(mylora, self).__init__(verbose)
self.set_mode(MODE.SLEEP)
self.set_dio_mapping([0] * 6)
def on_rx_done(self):
BOARD.led_on()
#print("\nRxDone")
self.clear_irq_flags(RxDone=1)
payload = self.read_payload(nocheck=True )# Receive INF
print ("Receive: ")
mens=bytes(payload).decode("utf-8",'ignore')
mens=mens[2:-1] #to discard \x00\x00 and \x00 at the end
print(mens)
BOARD.led_off()
if mens=="INF":
print("Received data request INF")
time.sleep(2)
print ("Send mens: DATA RASPBERRY PI")
self.write_payload([255, 255, 0, 0, 68, 65, 84, 65, 32, 82, 65, 83, 80, 66, 69, 82, 82, 89, 32, 80, 73, 0]) # Send DATA RASPBERRY PI
self.set_mode(MODE.TX)
time.sleep(2)
self.reset_ptr_rx()
self.set_mode(MODE.RXCONT)
def on_tx_done(self):
print("\nTxDone")
print(self.get_irq_flags())
def on_cad_done(self):
print("\non_CadDone")
print(self.get_irq_flags())
def on_rx_timeout(self):
print("\non_RxTimeout")
print(self.get_irq_flags())
def on_valid_header(self):
print("\non_ValidHeader")
print(self.get_irq_flags())
def on_payload_crc_error(self):
print("\non_PayloadCrcError")
print(self.get_irq_flags())
def on_fhss_change_channel(self):
print("\non_FhssChangeChannel")
print(self.get_irq_flags())
def start(self):
while True:
self.reset_ptr_rx()
self.set_mode(MODE.RXCONT) # Receiver mode
while True:
pass;
def send(data):
data = [255, 255, 0, 0] + data + [0]
self.write_payload([data]) # Send DATA
time.sleep(2)
self.reset_ptr_rx()
self.set_mode(MODE.RXCONT)
lora = mylora(verbose=False)
#args = parser.parse_args(lora) # configs in LoRaArgumentParser.py
# Slow+long range Bw = 125 kHz, Cr = 4/8, Sf = 4096chips/symbol, CRC on. 13 dBm
lora.set_pa_config(pa_select=1, max_power=21, output_power=15)
lora.set_bw(BW.BW125)
lora.set_coding_rate(CODING_RATE.CR4_8)
lora.set_spreading_factor(12)
lora.set_rx_crc(True)
#lora.set_lna_gain(GAIN.G1)
#lora.set_implicit_header_mode(False)
lora.set_low_data_rate_optim(True)
# Medium Range Defaults after init are 434.0MHz, Bw = 125 kHz, Cr = 4/5, Sf = 128chips/symbol, CRC on 13 dBm
#lora.set_pa_config(pa_select=1)
assert(lora.get_agc_auto_on() == 1)
try:
print("START")
lora.start()
except KeyboardInterrupt:
sys.stdout.flush()
print("Exit")
sys.stderr.write("KeyboardInterrupt\n")
finally:
sys.stdout.flush()
print("Exit")
lora.set_mode(MODE.SLEEP)
BOARD.teardown()
"""
signal.signal(signal.SIGINT, signal_handler)
def process_events(log=False,lorasend=False):
global cps
global cpm
nowtime = time.time()
prune_before_time = nowtime - start_time_epoch - 60.0
for pulse in list(pulse_events):
if(pulse < prune_before_time):
pulse_events.popleft()
cps = len(pulse_events)/60.0
cpm = len(pulse_events)
print(cpm)
if not(int(time.time()) % 60) and prune_before_time > 0:
# log last minute cpm to db. wait at least 60 sec from start
# to get steady state data.
print("last minute cpm:" + str(cpm))
if(lorasend):
epoch_ms = int(time.time()*1000.0)
buffer_data = struct.pack("<Q", epoch_ms) # pack epoch_ms into byte array little endian
buffer_data += struct.pack("<L", cpm) # append cpm into byte array little endian
#lora.send(buffer_data)
print(buffer_data)
if(log):
sql = "INSERT INTO data_cpm (timestamp_utc,count_per_minute,coordinates) VALUES(%s, %s, ST_PointFromWKB(ST_AsBinary(POINT(" + str(lat) + "," + str(long) + ")), 4326))";
val = (str(datetime.datetime.now()),str(cpm))
cur.execute(sql, val)
conn.commit()
processing_delay = time.time() - nowtime
print(processing_delay)
proc_delay_mul = int(processing_delay)
time.sleep(1 + proc_delay_mul - processing_delay)
# ensure isochronous sampling, add n skip seconds in case of the block taking a really long time...
# asynchronous mariadb cursor should prevent this occurence.
if(prune_before_time > 0): # 60 sec have elapsed
return cpm
else:
return -1
#while process_pulses:
# process_events()
#signal.pause()
Leave a Reply