Nuclear & radiation physicist: MQL recipes

Recipes

Activity of a check source today
Measure a half-life from counts
Mo-99 / Tc-99m generator: when to elute
Lead shielding for a Co-60 source
Buckling and peaking factor of a cylindrical core
Binding energy of the deuteron from the masses
Is the background really Poisson?
Neutron activation of a gold foil
Cosmic muons: time dilation in numbers
HPGe gamma peak: resolution, ROI and detection limit

executed on MentDB

Activity of a check source today

A Cs-137 check source was certified at 370 kBq on 2010-03-15. Compute its activity on 2026-11-01 (half-life 30.08 years), in kBq, in uCi and as a fraction of the certified value.

Pipeline

01 · halflifeHalf-life in daysunit convert→02 · elapsedDays since the certificatedate nb_day→03 · decayA = A0 2^(-t/T)math calc

Script

admin

#1. Cs-137 check source: 370 kBq on its calibration date, half-life 30.08 years;
-> "[A0]" 370;
-> "[T]" (unit convert 30.08 "yr" "d");
#2. Days elapsed since the certificate;
-> "[days]" (date nb_day "2010-03-15" "2026-11-01");
#3. Decay law A = A0 * 2^-t/T, 1 Ci = 3.7e10 Bq;
-> "[A]" (math calc "[A0]*2^(-[days]/[T])");
concat [days] " days elapsed, activity = " (math round [A] 1) " kBq = " (math round (math calc "[A]*1000/3.7e10*1e6") 3) " uCi, " (math round (math calc "[A]/[A0]*100") 1) " % of the certified value"

Result

mentdb

6075 days elapsed, activity = 252.2 kBq = 6.816 uCi, 68.2 % of the certified value

Go further

Use date sysdate to get the activity on the day you run it
Find when the source falls under 100 kBq with num solve

executed on MentDB

Measure a half-life from counts

Background-corrected counts of a Ba-137m sample are taken every minute. Fit N = N0 2^(-t/T) weighted by the Poisson errors sqrt(N), and report the half-life with its uncertainty and the reduced chi-square.

Pipeline

01 · fitWeighted fit of N0 and Tstat fit_model→02 · readParameters and errorsjson select→03 · reportT1/2 with its uncertaintymath round_uncertainty

Script

admin

#1. Ba-137m counts per 30 s, background subtracted, Poisson errors sqrt N;
-> "[t]" "[0, 1, 2, 3, 4, 5, 6]";
-> "[N]" "[5032, 3790, 2931, 2198, 1702, 1270, 991]";
-> "[s]" "[71, 62, 54, 47, 41, 36, 31]";
#2. Weighted fit of N = N0 * 2^-t/T;
json load "rx_hl" (stat fit_model "N0*2^(-t/T)" "t" [t] [N] "{\"N0\": 5000, \"T\": 2.5}" [s]);
#3. Half-life with its uncertainty and the goodness of fit;
concat "T1/2 = " (math round_uncertainty (json select "rx_hl" /parameters/T) (json select "rx_hl" /errors/T)) " min (tabulated 2.552 min), chi2/dof = " (math round (json select "rx_hl" /chi2_red) 2)

Result

mentdb

T1/2 = 2.55 ± 0.04 min (tabulated 2.552 min), chi2/dof = 0.21

Go further

Check the counting statistics with stat dist_cdf "poisson"
Compare with the linear fit of ln(N) using stat fit "exp"

executed on MentDB

Mo-99 / Tc-99m generator: when to elute

After an elution, Tc-99m grows again from Mo-99 (half-lives 65.94 h and 6.0067 h, branching 87.5 %). Maximise the Bateman activity of the daughter to find the best time for the next elution, and check it against t_max = ln(ld/lp)/(ld - lp).

Pipeline

01 · lambdaDecay constants ln2 / Tmath calc→02 · batemanDaughter activity A_d(t)concat→03 · maximizeTime of maximumnum maximize

Script

admin

#1. Mo-99 / Tc-99m generator: half-lives 65.94 h and 6.0067 h, branching 87.5 %;
-> "[lp]" (math calc "ln(2)/65.94"); -> "[ld]" (math calc "ln(2)/6.0067");
#2. Bateman: daughter activity for 1 unit of Mo-99 activity at elution time;
-> "[Ad]" (concat "0.875*" [ld] "/(" [ld] "-" [lp] ")*(exp(-" [lp] "*t)-exp(-" [ld] "*t))");
#3. Best time to elute again;
json load "rx_gen" (num maximize [Ad] "t" 0 100);
concat "next elution after " (math round (json select "rx_gen" /x) 2) " h (formula ln ld/lp / ld-lp = " (math round (math calc "ln([ld]/[lp])/([ld]-[lp])") 2) " h), Tc-99m activity = " (math round (math calc (concat (json select "rx_gen" /value) "*100")) 1) " % of the Mo-99 activity at the previous elution"

Result

mentdb

next elution after 22.84 h (formula ln ld/lp / ld-lp = 22.84 h), Tc-99m activity = 68.8 % of the Mo-99 activity at the previous elution

Go further

Solve the full chain with num ode {"Mo": "-lp*Mo", "Tc": "0.875*lp*Mo-ld*Tc"}
Tabulate the yield for a 24 h elution schedule with num table

executed on MentDB

Lead shielding for a Co-60 source

A worker stands 2 m from a 400 MBq Co-60 source (dose-rate constant 0.36 uSv m2/(MBq h)). Compute the unshielded dose rate, the half-value layer of lead at 1.25 MeV, and the lead thickness that brings the dose rate down to 2.5 uSv/h (narrow beam, no build-up).

Pipeline

01 · doseInverse square lawmath calc→02 · muAttenuation coefficient and HVLmath calc→03 · thicknessSolve H0 exp(-mu x) = 2.5num solve

Script

admin

#1. Co-60 source 400 MBq, dose-rate constant 0.36 uSv m2 per MBq h, worker at 2 m;
-> "[H0]" (math calc "0.36*400/2^2");
#2. Lead at 1.25 MeV: mu = mu/rho * rho = 0.0588 * 11.35 per cm;
-> "[mu]" (math calc "0.0588*11.35");
#3. Lead thickness for 2.5 uSv/h, no build-up factor;
-> "[x]" (num solve (concat [H0] "*exp(-" [mu] "*x)-2.5") "x" 0 20);
concat "unshielded: " [H0] " uSv/h; HVL = " (math round (math calc "ln(2)/[mu]") 2) " cm; " (math round [x] 2) " cm of lead = " (math round (math calc "[x]*[mu]/ln(2)") 1) " HVL bring it to 2.5 uSv/h"

Result

mentdb

unshielded: 36 uSv/h; HVL = 1.04 cm; 4 cm of lead = 3.8 HVL bring it to 2.5 uSv/h

Go further

Add a linear build-up factor B = 1 + mu x and solve again with num solve
Find the distance giving the same dose rate without shield with num solve

executed on MentDB

Buckling and peaking factor of a cylindrical core

In a bare cylindrical reactor the flux is J0(2.405 r/R) cos(pi z/H). The script finds the first zero of J0, the geometric buckling, and the peak-to-average flux ratio by a double integral over the core, checked against the analytical value.

Pipeline

01 · zeroFirst zero of J0num solve→02 · bucklingGeometric bucklingmath calc→03 · integralFlux over the core volumenum integrate2→04 · peakingPeak to average ratiobessel_j(nu, x)

Script

admin

#1. Bare cylindrical core R = 1 m, H = 2 m;
-> "[R]" 1; -> "[H]" 2;
#2. First zero of J0 and geometric buckling;
-> "[x0]" (num solve "bessel_j(0, x)" "x" 2 3);
-> "[B2]" (math calc "([x0]/[R])^2+(pi/[H])^2");
#3. Flux shape J0 of x0 r/R * cos of pi z/H: volume integral over the core;
-> "[I]" (num integrate2 "2*pi*r*bessel_j(0, 2.404825557695773*r)*cos(pi*z/2)" "r" 0 1 "z" -1 1);
#4. Peak to average flux ratio, numerical and analytical;
-> "[ratio]" (math calc "pi*[R]^2*[H]/[I]");
-> "[theory]" (math calc "[x0]/(2*bessel_j(1, [x0]))*pi/2");
concat "{\"x0\": " [x0] ", \"buckling_m-2\": " (math round [B2] 4) ", \"peak_to_average\": " (math round [ratio] 4) ", \"analytical\": " (math round [theory] 4) "}"

Result

mentdb

{"x0": 2.404825557695773, "buckling_m-2": 8.2506, "peak_to_average": 3.6382, "analytical": 3.6382}

Go further

Minimize the core volume for a given buckling with num minimize: the optimum is H/D = 0.924.
A sphere uses sph_bessel_j(0, x), whose first zero is pi.

executed on MentDB

Binding energy of the deuteron from the masses

The mass defect of the deuteron, delta m = mp + mn - md, times c^2 gives its binding energy (2.2246 MeV, the energy of the gamma ray emitted in n + p capture). Also the binding energy per nucleon and the mass defect in atomic mass units.

Pipeline

01 · massesProton, neutron, deuteronconstant physics_m_deuteron→02 · c2c squaredconstant physics_c_square→03 · energyJoules to MeVunit convert→04 · amuDefect in uconstant physics_amu

Script

admin

-> "[mp]" (string regexp_extract (constant physics_m_proton) "[-0-9.]+E?[-0-9]*" 0);
-> "[mn]" (string regexp_extract (constant physics_m_neutron) "[-0-9.]+E?[-0-9]*" 0);
-> "[md]" (string regexp_extract (constant physics_m_deuteron) "[-0-9.]+E?[-0-9]*" 0);
-> "[c2]" (string regexp_extract (constant physics_c_square) "[-0-9.]+E?[-0-9]*" 0);
-> "[u]" (string regexp_extract (constant physics_amu) "[-0-9.]+E?[-0-9]*" 0);
#1. Mass defect in kg;
-> "[dm]" (math calc "[mp]+[mn]-[md]");
#2. Binding energy;
-> "[B]" (unit convert (math calc "[dm]*[c2]") "J" "MeV");
concat "delta m = " (math round_sig (math calc "[dm]*1E30") 5) "E-30 kg = " (math round_sig (math calc "[dm]/[u]") 5) " u; B = " (math round [B] 4) " MeV, B/A = " (math round (math calc "[B]/2") 4) " MeV; capture gamma wavelength " (math round (math calc (concat "1239.84198/" [B])) 2) " fm"

Result

mentdb

delta m = 3.9656E-30 kg = 0.0023882 u; B = 2.2246 MeV, B/A = 1.1123 MeV; capture gamma wavelength 557.34 fm

Go further

Do the same for helium-4 (28.3 MeV): why is it so tightly bound?
The CODATA 2022 masses have 11 to 12 digits: const info "mp" gives their uncertainty.

executed on MentDB

Is the background really Poisson?

Thirty 10-second background counts from a GM tube. Frequency table, mode, mean and population variance: for a Poisson process the dispersion index var/mean is close to 1. Then the Poisson probability of the modal count, lambda^k e^-lambda / Gamma(k + 1).

Pipeline

01 · tableFrequency tablestat frequency→02 · modeMost frequent countstat mode→03 · dispersionVariance over meanstat pvariance→04 · poissonPoisson probabilitymath gamma

Script

admin

-> "[c]" "[2, 7, 5, 3, 9, 4, 6, 1, 5, 8, 4, 6, 3, 7, 5, 2, 6, 4, 10, 5, 3, 6, 4, 7, 5, 8, 3, 5, 6, 2]";
#1. Centre and spread;
-> "[mean]" (math avg [c]);
-> "[var]" (stat pvariance [c]);
-> "[k]" (string regexp_extract (stat mode [c]) "[0-9]+" 0);
-> "[freq]" (stat frequency [c]);
-> "[obs]" (string regexp_extract [freq] (concat "\"value\": " [k] ", \"count\": [0-9]+, \"frequency\": ([0-9.]+)") 1);
#2. Poisson probability of the mode;
-> "[p]" (math calc (concat (math power [mean] [k]) "*exp(-[mean])/" (math gamma (math calc "[k]+1"))));
concat "mean = " [mean] ", variance = " (math round [var] 4) ", dispersion index = " (math round (math calc "[var]/[mean]") 3) "; mode = " [k] ", observed frequency " (math round [obs] 3) ", Poisson " (math round [p] 3) "; table " [freq]

Result

mentdb

mean = 5.033333333333333, variance = 4.6322, dispersion index = 0.92; mode = 5, observed frequency 0.2, Poisson 0.175; table [
  {"value": 1, "count": 1, "frequency": 0.03333333333333333, "cumulative_count": 1, "cumulative_frequency": 0.03333333333333333},
  {"value": 2, "count": 3, "frequency": 0.1, "cumulative_count": 4, "cumulative_frequency": 0.13333333333333333},
  {"value": 3, "count": 4, "frequency": 0.13333333333333333, "cumulative_count": 8, "cumulative_frequency": 0.26666666666666666},
  {"value": 4, "count": 4, "frequency": 0.13333333333333333, "cumulative_count": 12, "cumulative_frequency": 0.4},
  {"value": 5, "count": 6, "frequency": 0.2, "cumulative_count": 18, "cumulative_frequency": 0.6},
  {"value": 6, "count": 5, "frequency": 0.16666666666666666, "cumulative_count": 23, "cumulative_frequency": 0.7666666666666667},
  {"value": 7, "count": 3, "frequency": 0.1, "cumulative_count": 26, "cumulative_frequency": 0.8666666666666667},
  {"value": 8, "count": 2, "frequency": 0.06666666666666667, "cumulative_count": 28, "cumulative_frequency": 0.9333333333333333},
  {"value": 9, "count": 1, "frequency": 0.03333333333333333, "cumulative_count": 29, "cumulative_frequency": 0.9666666666666667},
  {"value": 10, "count": 1, "frequency": 0.03333333333333333, "cumulative_count": 30, "cumulative_frequency": 1}
]

Go further

A dispersion index far below 1 can betray dead-time losses; far above, an unstable source or electronics noise.
Test it formally: 29 x index follows a chi-square with 29 degrees of freedom, see stat dist_cdf.

executed on MentDB

Neutron activation of a gold foil

A 0.1 g gold foil (Au-197, sigma = 98.65 b) sits 2 h in a thermal flux of 1E12 n/cm^2.s, then cools 1 h. Au-198 has T1/2 = 2.695 d. A = N sigma phi (1 - e^(-lambda t_irr)) e^(-lambda t_d), with 1 - e^(-x) computed by expm1.

Pipeline

01 · atomsAtoms with Avogadroconstant physics_n→02 · lambdaDecay constantunit convert→03 · growthSaturation factormath expm1→04 · decayCoolingmath exp

Script

admin

-> "[NA]" (string regexp_extract (constant physics_n) "[-0-9.]+E?[-0-9]*" 0);
#1. Atoms of Au-197 and reaction rate;
-> "[Nt]" (math calc "0.1/196.967*[NA]");
-> "[R]" (math calc "[Nt]*98.65E-24*1E12");
#2. Decay constant of Au-198;
-> "[l]" (math calc (concat "ln(2)/" (unit convert 2.695 "d" "s")));
#3. Build-up during 2 h then 1 h of decay;
-> "[S]" (math calc (concat "-" (math expm1 (math calc "-[l]*7200"))));
-> "[A]" (math calc (concat "[R]*[S]*" (math exp (math calc "-[l]*3600"))));
concat "N = " (math round_sig (math calc "[Nt]/1E20") 4) "E20 atoms, reaction rate " (math round_sig (math calc "[R]/1E10") 4) "E10 1/s; saturation factor " (math round_sig [S] 4) "; activity after cooling " (math round (math calc "[A]/1E6") 2) " MBq = " (math round (math calc "[A]/3.7E7") 2) " mCi; saturation activity " (math round (math calc "[R]/1E9") 3) " GBq"

Result

mentdb

N = 3.057E20 atoms, reaction rate 3.016E10 1/s; saturation factor 0.02121; activity after cooling 632.76 MBq = 17.1 mCi; saturation activity 30.162 GBq

Go further

Measure the activity and invert the formula to get the flux: this is how flux maps are made.
Cadmium-covered foils separate the epithermal part (cadmium ratio).

executed on MentDB

Cosmic muons: time dilation in numbers

A muon created 15 km up with a total energy of 4 GeV. Rest energy from the muon mass, gamma = E/mc^2, rapidity acosh(gamma), beta = tanh(rapidity), mean decay length beta gamma c tau, and the fraction reaching the ground with and without relativity.

Pipeline

01 · massMuon rest energyconstant physics_m_muon→02 · rapidityRapidity from gammamath acosh→03 · betaVelocitymath tanh→04 · survivalExponential decaymath exp

Script

admin

-> "[c2]" (string regexp_extract (constant physics_c_square) "[-0-9.]+E?[-0-9]*" 0);
-> "[mc2]" (unit convert (math calc (concat (string regexp_extract (constant physics_m_muon) "[-0-9.]+E?[-0-9]*" 0) "*[c2]")) "J" "MeV");
#1. Kinematics;
-> "[g]" (math calc "4000/[mc2]");
-> "[y]" (math acosh [g]);
-> "[b]" (math tanh [y]);
#2. Mean decay length with tau = 2.197 us;
-> "[L]" (math calc "[b]*[g]*299792458*2.197E-6");
concat "mc2 = " (math round [mc2] 3) " MeV, gamma = " (math round [g] 2) ", rapidity = " (math round [y] 4) ", beta = " (math round [b] 6) "; decay length " (math round (math calc "[L]/1000") 2) " km; survival over 15 km: " (math round (math calc (concat "100*" (math exp (math calc "-15000/[L]")))) 1) " % with relativity, " (math round_sig (math calc (concat "100*" (math exp (math calc "-15000/(299792458*2.197E-6)")))) 2) " % without"

Result

mentdb

mc2 = 105.658 MeV, gamma = 37.86, rapidity = 4.3268, beta = 0.999651; decay length 24.93 km; survival over 15 km: 54.8 % with relativity, 0.000000013 % without

Go further

The real muon spectrum peaks near 4 GeV at sea level: integrate over it with num integrate.
In the muon frame the atmosphere is Lorentz-contracted by the same gamma.

executed on MentDB

HPGe gamma peak: resolution, ROI and detection limit

The Co-60 1332.5 keV peak has a FWHM of 1.8 keV. sigma = FWHM / (2 sqrt(2 ln 2)); fraction of the peak inside +-1 FWHM with erf; ROI half-width that keeps 99.7 % with erf_inv; Currie detection limit L_D = 2.71 + 4.65 sqrt(B) for 400 background counts in the ROI.

Pipeline

01 · sigmaFWHM to sigmamath log→02 · fractionCounts inside the ROImath erf→03 · roiROI for 99.7 %math erf_inv→04 · currieDetection limitmath sqrt

Script

admin

-> "[fw]" 1.8;
#1. Gaussian width;
-> "[s]" (math calc (concat "[fw]/(2*sqrt(2*" (math log 2) "))"));
#2. Fraction within +- 1 FWHM;
-> "[f]" (math erf (math calc "[fw]/([s]*sqrt(2))"));
#3. Half-width for 99.7 percent;
-> "[hw]" (math calc (concat (math erf_inv 0.997) "*sqrt(2)*[s]"));
#4. Currie limit with 400 background counts;
-> "[LD]" (math calc (concat "2.71+4.65*" (math sqrt 400)));
concat "sigma = " (math round [s] 4) " keV, resolution " (math round (math calc "100*[fw]/1332.5") 3) " %; +-1 FWHM keeps " (math round (math calc "100*[f]") 2) " %; 99.7 % needs +-" (math round [hw] 3) " keV = " (math round (math calc "[hw]/[s]") 3) " sigma; L_D = " (math round [LD] 1) " counts"

Result

mentdb

sigma = 0.7644 keV, resolution 0.135 %; +-1 FWHM keeps 98.15 %; 99.7 % needs +-2.269 keV = 2.968 sigma; L_D = 95.7 counts

Go further

A wider ROI keeps more signal but adds background: the optimum for a weak peak is about +-1.2 FWHM.
Fit the peak with stat fit_model and a Gaussian on a linear background.