Medical physicist: MQL recipes

Recipes

I-131 therapy: decay correction and MIRD thyroid dose
Beam quality (HVL) and dual-energy separation
MTF of a detector from its line spread function
Hypofractionation with the linear-quadratic model and TCP
Linac bunker: primary barrier thickness (NCRP 147)
Triple-isotope SPECT: crosstalk correction by LU
HDR brachytherapy: TG-43 line-source dose rate
MRI: Larmor frequency, RF wavelength and STIR null point
Well counter constancy: chi-square and z-test
Proton therapy: energy layers with the Bragg-Kleeman rule

executed on MentDB

I-131 therapy: decay correction and MIRD thyroid dose

A 400 MBq I-131 capsule calibrated two days before the treatment of a 25 g hyperthyroid gland with 45 % uptake. The script corrects the activity for decay, computes the effective half-life, integrates the decay curve exactly with the CAS, and gives the cumulated activity, a beta S value from the mean beta energy, the absorbed dose and its uncertainty from the uptake and the biological half-life.

Pipeline

01 · decayA = A0 2^(-t/T)math calc→02 · TeEffective half-lifebigdec round→03 · integralIntegral of exp(-ln2 t / Te) from 0 to infinitycas integrate · math eval→04 · S value0.1916 MeV / 25 g in Gyunit calc · unit convert · unit prefix→05 · doseD = A~ x S and its uncertaintybigdec precision · num uncertainty

Script

admin

#1. I-131 capsule for a hyperthyroid patient, calibrated at 400 MBq two days before;
-> "[Tp]" 8.02; -> "[A0]" (math calc "400*2^(-2/[Tp])");
#2. Effective half-life with a biological half-life of 20 days;
-> "[Te]" (bigdec round (math calc "[Tp]*20/([Tp]+20)") 3);
#3. Integral of the decay curve from 0 to infinity: exact with the CAS, then in seconds;
-> "[ex]" (cas integrate (concat "exp(-ln(2)*t/" [Te] ")") "t" 0 "Infinity");
-> "[Ac]" (math calc (concat "0.45*[A0]*" (math eval [ex]) "*" (unit convert 1 "d" "s")));
#4. MIRD S value: mean beta energy 0.1916 MeV absorbed in a 25 g gland, in mGy per MBq.s;
-> "[S]" (math calc (concat (unit calc "0.1916 MeV / 25 g" "Gy") "*1e6*1e3"));
-> "[D]" (math calc "[Ac]*[S]/1000");
#5. Uncertainty: uptake 0.45 +/- 0.05, biological half-life 20 +/- 5 days;
-> "[u]" (json load select "rx_u" (num uncertainty "f*A*(T*Tb/(T+Tb))/ln(2)*86400*S/1000" (concat "{\"f\": [0.45, 0.05], \"A\": " [A0] ", \"T\": " [Tp] ", \"Tb\": [20, 5], \"S\": " [S] "}")) "/text");
concat "administered " (bigdec precision [A0] 4) " MBq, Te = " [Te] " d, integral " [ex] " d, cumulated activity " (bigdec precision (unit prefix [Ac] M G) 3) " GBq.s, S = " (bigdec precision [S] 3) " mGy/(MBq.s), thyroid dose = " (bigdec precision [D] 3) " Gy, with uncertainties " [u] " Gy"

Result

mentdb

administered 336.5 MBq, Te = 5.724 d, integral 1431/250*1/ln(2) d, cumulated activity 108000 GBq.s, S = 0.00123 mGy/(MBq.s), thyroid dose = 133 Gy, with uncertainties 133 ± 18 Gy

Go further

Measure the biological half-life from uptake points at 4, 24 and 96 h with stat fit_model
Solve the activity that gives a target dose of 150 Gy with num solve
Add the gamma contribution with a second S value

executed on MentDB

Beam quality (HVL) and dual-energy separation

From the transmission of an X-ray beam through aluminium sheets, fit the effective attenuation coefficient, interpolate the first and second half-value layers and the homogeneity coefficient. Then separate soft tissue and bone thicknesses from the log attenuation measured at two energies (2 x 2 linear system) and check how well the system is conditioned.

Pipeline

01 · fitT = exp(-mu x)stat fit_model→02 · HVLx at T = 0.5 and 0.25num interpolate→03 · separateSolve mu . t = ln(I0/I)matrix solve→04 · checkCondition number and determinantmatrix cond · matrix det

Script

admin

#1. Transmission of the beam through aluminium sheets, thickness in mm;
-> "[x]" "[0, 1, 2, 3, 4, 6, 8]";
-> "[T]" "[1, 0.781, 0.618, 0.497, 0.405, 0.276, 0.193]";
#2. Effective attenuation coefficient and first half-value layer;
-> "[mu]" (json load select "rx_mu" (stat fit_model "exp(-mu*x)" "x" [x] [T] "{\"mu\": 0.2}") "/parameters/mu");
-> "[hvl1]" (num interpolate [T] [x] 0.5 "cubic");
-> "[hvl2]" (math calc (concat (num interpolate [T] [x] 0.25 "cubic") "-[hvl1]"));
#3. Dual-energy separation: ln I0/I = mu_soft * t_soft + mu_bone * t_bone at 2 energies, in cm-1 and cm;
-> "[M]" "[[0.227, 0.573], [0.184, 0.302]]";
-> "[t]" (matrix solve [M] "[4.79, 3.42]");
concat "mu_eff = " (math round [mu] 4) " /mm, HVL1 = " (math round [hvl1] 2) " mm Al, HVL2 = " (math round [hvl2] 2) " mm Al, homogeneity coefficient " (math round (/ [hvl1] [hvl2]) 3) "; soft tissue and bone thickness " [t] " cm, condition number " (math round (matrix cond [M]) 1) ", det " (math round (matrix det [M]) 4)

Result

mentdb

mu_eff = 0.2248 /mm, HVL1 = 2.97 mm Al, HVL2 = 3.58 mm Al, homogeneity coefficient 0.83; soft tissue and bone thickness [13.912902001193123,2.847768317153859] cm, condition number 13.6, det -0.0369

Go further

Propagate a 1 % noise on the measured log attenuations with num uncertainty or matrix inverse
Plot the transmission data with plot points
Compare the HVL to the minimum required for the tube voltage

executed on MentDB

MTF of a detector from its line spread function

The line spread function measured every 0.1 mm (Gaussian, sigma 0.25 mm, with a small ripple) is transformed with an FFT; the modulus normalized at f = 0 is the MTF. The script interpolates the MTF50, compares it with the Gaussian theory and with the aperture of a 2 mm pinhole, MTF = 2 J1(pi d f) / (pi d f), and draws both curves.

Pipeline

01 · LSF64 samples in a loopmath calc→02 · FFTFourier transform and modulusnum fft · complex abs→03 · MTF50Interpolate f at MTF = 0.5num interpolate→04 · pinholeRoot of 2 J1(x)/x = 0.5num solve · bessel_j(nu, x)→05 · plotMeasured and theoretical MTFplot chart

Script

admin

#1. Line spread function: 64 samples every 0.1 mm, Gaussian sigma 0.25 mm plus a small ripple;
-> "[lsf]" "["; -> "[sep]" "";
for (-> "[i]" 0) (< [i] 64) (++ "[i]") {
  -> "[lsf]" (concat [lsf] [sep] (math calc "exp(-([i]*0.1-3.2)^2/(2*0.25^2))+0.002*sin(1.7*[i])"));
  -> "[sep]" ",";
};
-> "[lsf]" (concat [lsf] "]");
#2. MTF = modulus of the Fourier transform normalized at f = 0, pinhole aperture 2 J1 of pi d f over pi d f;
json load "rx_ft" (num fft [lsf]);
-> "[m0]" (complex abs (json select "rx_ft" "/[0]"));
-> "[f]" "[0"; -> "[mtf]" "[1"; -> "[ph]" "[1";
for (-> "[k]" 1) (< [k] 17) (++ "[k]") {
  -> "[f]" (concat [f] "," (math calc "[k]/6.4"));
  -> "[mtf]" (concat [mtf] "," (math calc (concat (complex abs (json select "rx_ft" (concat "/[" [k] "]"))) "/[m0]")));
  -> "[ph]" (concat [ph] "," (math calc "2*bessel_j(1, pi*2*[k]/6.4)/(pi*2*[k]/6.4)"));
};
-> "[f]" (concat [f] "]"); -> "[mtf]" (concat [mtf] "]"); -> "[ph]" (concat [ph] "]");
#3. MTF50 of the detector and of a 2 mm pinhole;
-> "[m50]" (num interpolate [mtf] [f] 0.5 "cubic");
-> "[p50]" (num solve "2*bessel_j(1, pi*2*x)/(pi*2*x)-0.5" "x" 0.1 0.6);
plot chart (concat "{\"title\": \"MTF50: detector " (math round [m50] 3) " lp/mm, Gaussian theory " (math round (math calc "sqrt(ln(2)/2)/(pi*0.25)") 3) ", pinhole 2 mm " (math round [p50] 3) " lp/mm\", \"x_label\": \"spatial frequency (lp/mm)\", \"y_label\": \"MTF\", \"series\": [{\"name\": \"measured MTF\", \"type\": \"points\", \"x\": " [f] ", \"y\": " [mtf] "}, {\"name\": \"pinhole aperture\", \"x\": " [f] ", \"y\": " [ph] "}]}")

Result

mentdb

The chart is shown in the Plot tab of the editor.

Go further

Apply dsp window "hann" to the tails of a noisy LSF before the FFT
Get the noise power spectrum of a flat image row with dsp welch
Combine detector and aperture: the system MTF is the product

executed on MentDB

Hypofractionation with the linear-quadratic model and TCP

Convert a conventional prostate treatment (39 x 2 Gy) into a 5-fraction scheme with the same tumour BED (alpha/beta = 1.5 Gy), compare the late rectal EQD2 (alpha/beta = 3 Gy), and compute the Poisson tumour control probability, also for a de-escalated 5 x 6 Gy.

Pipeline

01 · BEDn.d(1 + d/(alpha/beta))math calc→02 · iso-effectDose per fraction in 5 fractionsnum solve→03 · EQD2Rectum, alpha/beta = 3math calc→04 · TCPP(0 surviving clonogens)stat dist_pdf

Script

admin

#1. Prostate: conventional 39 x 2 Gy against a hypofractionated scheme in 5 fractions, alpha/beta = 1.5 Gy tumour, 3 Gy rectum;
-> "[ab]" 1.5; -> "[abr]" 3;
-> "[BED]" (math calc "39*2*(1+2/[ab])");
#2. Dose per fraction d in 5 fractions giving the same tumour BED, linear-quadratic model;
-> "[d]" (num solve (concat "5*x*(1+x/" [ab] ")-" [BED]) "x" 2 15);
#3. Equivalent doses in 2 Gy fractions, tumour and late rectal effects;
-> "[eqd2r1]" (math calc "39*2*(2+[abr])/(2+[abr])");
-> "[eqd2r2]" (math calc "5*[d]*([d]+[abr])/(2+[abr])");
#4. Poisson tumour control probability: 1e7 clonogens, alpha = 0.10 per Gy, SF = exp of -alpha.BED, TCP = P of zero survivors;
-> "[tcp]" (stat dist_pdf "poisson" (concat "[" (math calc "1e7*exp(-0.10*[BED])") "]") 0);
#Same with a de-escalated 5 x 6 Gy;
-> "[tcp6]" (stat dist_pdf "poisson" (concat "[" (math calc "1e7*exp(-0.10*5*6*(1+6/[ab]))") "]") 0);
concat "BED = " (math round [BED] 1) " Gy, iso-effective dose = 5 x " (math round [d] 2) " Gy = " (math round (* 5 [d]) 1) " Gy; rectum EQD2 " (math round [eqd2r1] 1) " Gy conventional vs " (math round [eqd2r2] 1) " Gy hypofractionated; TCP = " (math round [tcp] 3) ", with 5 x 6 Gy only " (math round [tcp6] 3)

Result

mentdb

BED = 182 Gy, iso-effective dose = 5 x 6.68 Gy = 33.4 Gy; rectum EQD2 78 Gy conventional vs 64.6 Gy hypofractionated; TCP = 0.883, with 5 x 6 Gy only 0.047

Go further

Add the repopulation term (T - Tk).ln2/(alpha.Tpot) for a prolonged treatment
Draw TCP against the dose per fraction with plot function

executed on MentDB

Linac bunker: primary barrier thickness (NCRP 147)

Design the primary barrier of a 6 MV linac: distance from the target to the protected point, transmission B = P.d^2/(W.U.T) for a controlled area, number of tenth-value layers and thickness of ordinary concrete.

Pipeline

01 · distance|target - point|vector norm→02 · transmissionB = P.d^2/(W.U.T)math calc→03 · TVLn = log10(1/B)math log10→04 · thicknessTVL1 + (n-1).TVLeunit convert

Script

admin

#1. Primary barrier of a 6 MV linac bunker: workload 450 Gy per week at 1 m, use factor 0.25, controlled area 0.1 mSv per week, occupancy 1;
-> "[P]" 0.0001; -> "[W]" 450; -> "[U]" 0.25; -> "[T]" 1;
#2. Distance from the target to the protected point 0.3 m beyond the wall, coordinates in m;
-> "[dist]" (vector norm (vector sub "[0, 5.9, 1.0]" "[0, -0.4, 0.2]"));
#3. NCRP 147: transmission B = P.d^2 over W.U.T, number of tenth-value layers;
-> "[B]" (math calc "[P]*[dist]^2/([W]*[U]*[T])");
-> "[n]" (math log10 (/ 1 [B]));
#4. Ordinary concrete 2.35 g/cm3: first TVL 37 cm, equilibrium TVL 33 cm;
-> "[tc]" (math calc "37+([n]-1)*33");
concat "d = " (math round [dist] 3) " m, B = " (math round_sig [B] 3) ", n = " (math round [n] 3) " TVL; concrete " (math round [tc] 1) " cm = " (math round (unit convert [tc] "cm" "ft") 2) " ft, built as " (* 5 (math ceiling (/ [tc] 5))) " cm"

Result

mentdb

d = 6.351 m, B = 0.0000358, n = 4.446 TVL; concrete 150.7 cm = 4.94 ft, built as 155 cm

Go further

Do the secondary barrier: leakage at 0.1 % and scatter with the patient scatter fraction
Check the instantaneous dose rate at the point and the time-averaged dose rate

executed on MentDB

Triple-isotope SPECT: crosstalk correction by LU

In a Tc-99m / I-123 / In-111 acquisition, each energy window also counts the photons of the other isotopes. With the window sensitivities measured on phantoms, an LU decomposition separates the three activities (forward and back substitution).

Pipeline

01 · decomposeP.A = L.Umatrix lu→02 · permuteP.bmatrix mul→03 · substituteL.y = P.b then U.x = ymatrix solve→04 · determinantProduct of the U diagonalmatrix get

Script

admin

#1. Triple-isotope SPECT: counts per MBq of Tc-99m, I-123, In-111 rows in the 140, 159 and 245 keV windows;
-> "[a]" "[[820, 260, 40], [150, 710, 95], [10, 30, 640]]";
-> "[b]" "[21400, 17350, 6990]";
#2. LU decomposition with partial pivoting P.A = L.U;
json load "rx_p4d_lu" (matrix lu [a]);
-> "[l]" (json select "rx_p4d_lu" "/L");
-> "[u]" (json select "rx_p4d_lu" "/U");
-> "[p]" (json select "rx_p4d_lu" "/P");
#3. Forward then back substitution: L.y = P.b, U.x = y;
-> "[y]" (matrix solve [l] (matrix mul [p] [b]));
-> "[x]" (matrix solve [u] [y]);
#4. Crosstalk-free activities and the gain against the raw photopeak estimate;
concat "activities Tc, I, In = " [x] " MBq, |det A| = product of the U diagonal = " (math abs (math round (* (matrix get [u] 1 1) (* (matrix get [u] 2 2) (matrix get [u] 3 3))) 0)) ", raw Tc estimate " (math round (/ 21400 820) 2) " MBq"

Result

mentdb

activities Tc, I, In = [19.600545948230444,18.994392885883503,9.72525430303311] MBq, |det A| = product of the U diagonal = 3.45454E8, raw Tc estimate 26.1 MBq

Go further

Check the conditioning of the sensitivity matrix with matrix cond: close windows amplify the noise.
Propagate the Poisson noise of the counts with (A^-1)^2 . b (matrix inverse).

executed on MentDB

HDR brachytherapy: TG-43 line-source dose rate

For an Ir-192 HDR source, the script gets the polar coordinates of a calculation point, the angle subtended by the active length, the TG-43 line-source geometry factor normalised at the reference point, and the dose rate.

Pipeline

01 · radiusr = hypot(x, z)math hypot→02 · angletheta and betamath atan2→03 · referenceG_L(1 cm, 90 deg)math atan→04 · geometrybeta / (L r sin theta)math sin→05 · doseSk.Lambda.G.g.Fmath calc

Script

admin

#1. HDR Ir-192 line source of 3.5 mm along z, air kerma strength 40 700 U, point at x = 1.5 cm, z = 2 cm;
-> "[sk]" 40700; -> "[lam]" 1.108; -> "[len]" 0.35; -> "[x]" 1.5; -> "[z]" 2;
#2. Polar coordinates of the point: r with hypot, theta from the source axis with atan2;
-> "[r]" (math hypot [x] [z]);
-> "[th]" (math atan2 [x] [z]);
#3. Angle beta subtended by the source, from both ends of the line;
-> "[beta]" (- (math atan2 [x] (- [z] (/ [len] 2))) (math atan2 [x] (+ [z] (/ [len] 2))));
#4. TG-43 line-source geometry factors, here and at the reference point r0 = 1 cm, theta0 = 90 deg;
-> "[g]" (/ [beta] (* [len] (* [r] (math sin [th]))));
-> "[g0]" (/ (* 2 (math atan (/ [len] 2))) [len]);
#5. Dose rate with g of r = 0.99 and F of r,theta = 0.97 from the consensus tables;
-> "[d]" (* [sk] (* [lam] (* (/ [g] [g0]) (* 0.99 0.97))));
concat "r = " [r] " cm, theta = " (math round (math calc "[th]*180/pi") 2) " deg, G_L ratio = " (math round (/ [g] [g0]) 5) " vs 1/r2 = " (math round (/ 1 (* [r] [r])) 5) ", dose rate = " (math round (/ [d] 3600) 2) " cGy/s"

Result

mentdb

r = 2.5 cm, theta = 36.87 deg, G_L ratio = 0.16203 vs 1/r2 = 0.16, dose rate = 1.95 cGy/s

Go further

Interpolate g(r) and F(r, theta) in the consensus tables with num interpolate2.
Sum the dwell positions of a catheter with a for loop for a full plan point.

executed on MentDB

MRI: Larmor frequency, RF wavelength and STIR null point

From the CODATA proton magnetic moment the script computes the proton gyromagnetic ratio, the Larmor frequencies at 1.5 T and 3 T, the RF wavelength in tissue (dielectric effects at 3 T), the STIR inversion time that nulls fat and the white matter signal left.

Pipeline

01 · gamma2.mu_p/hmath calc→02 · referencePlanck constantconstant physics_planck→03 · wavelengthc / (f.sqrt(eps_r))math calc→04 · STIRTI = T1.ln 2math log→05 · signal|1 - 2e^(-TI/T1)|math exp

Script

admin

#1. Proton Larmor frequency f = 2.mu_p.B/h from the CODATA constants;
-> "[gbar]" (math calc "2*_mup/_h/1e6");
#2. At 1.5 T and 3 T, and the RF wavelength in tissue with a relative permittivity of 50;
-> "[f15]" (* [gbar] 1.5); -> "[f3]" (* [gbar] 3);
-> "[lam3]" (math calc "_c/([f3]*1e6*sqrt(50))*100");
#3. STIR fat suppression: null point TI = T1.ln 2, fat T1 = 260 ms at 1.5 T and 380 ms at 3 T;
-> "[ti15]" (* 260 (math log 2)); -> "[ti3]" (* 380 (math log 2));
#4. Signal left for white matter T1 = 1080 ms at 3 T with this TI: abs of 1 - 2.e^-TI/T1;
-> "[s]" (math abs (- 1 (* 2 (math exp (/ (- 0 [ti3]) 1080)))));
concat "gamma/2pi = " (math round [gbar] 4) " MHz/T (h = " (constant physics_planck) "), f = " (math round [f15] 2) " MHz and " (math round [f3] 2) " MHz, RF wavelength in tissue at 3 T = " (math round [lam3] 1) " cm, STIR TI = " (math round [ti15] 0) " ms / " (math round [ti3] 0) " ms, white matter signal " (math round (* 100 [s]) 1) " % of M0"

Result

mentdb

gamma/2pi = 42.5775 MHz/T (h = 6.62607015E-34 J·s), f = 63.87 MHz and 127.73 MHz, RF wavelength in tissue at 3 T = 33.2 cm, STIR TI = 180 ms / 263 ms, white matter signal 56.7 % of M0

Go further

Compute the SAR scaling with B0^2 between 1.5 T and 3 T.
Find the TI that nulls the CSF in FLAIR (T1 = 4000 ms) and the signal of the lesions.

executed on MentDB

Well counter constancy: chi-square and z-test

The daily quality control of a gamma well counter repeats 10 counts of the Cs-137 check source. The script decay-corrects the expected count, runs the Poisson chi-square test (accepted between 0.1 and 0.9) and a z-test of accuracy with sigma = sqrt(N).

Pipeline

01 · decayA0.exp(-lambda t)math exp→02 · PoissonPopulation variancestat pvariance→03 · chi2P(chi2(9))stat dist_cdf→04 · accuracyz-teststat ztest

Script

admin

#1. Well counter constancy: 10 counts of 60 s on the Cs-137 check source, background subtracted;
-> "[c]" "[48210, 47880, 48450, 48120, 47990, 48430, 48300, 47950, 48540, 48070]";
#2. Expected count: 50 500 at acceptance, decayed over 2.1 years, Cs-137 half-life 30.08 years;
-> "[exp]" (* 50500 (math exp (* (- 0 (/ (math log 2) 30.08)) 2.1)));
#3. Poisson check: the population variance should be close to the mean, chi2 = sum of x-mean2 / mean;
-> "[m]" (stat mean [c]);
-> "[chi2]" (/ (* 10 (stat pvariance [c])) [m]);
-> "[pchi]" (stat dist_cdf "chi2" "[9]" [chi2]);
#4. Accuracy: z-test of the mean against the expected count with sigma = sqrt of the count;
-> "[zt]" (stat ztest [c] [exp] (math sqrt [exp]));
concat "expected = " (math round [exp] 0) ", mean = " [m] " (" (math round (* 100 (- (/ [m] [exp]) 1)) 2) " %), chi2 = " (math round [chi2] 2) " with P chi2 <= = " (math round [pchi] 3) " (accept 0.1-0.9), z-test " [zt]

Result

mentdb

expected = 48114, mean = 48194 (0.17 %), chi2 = 9.82 with P chi2 <= = 0.634 (accept 0.1-0.9), z-test {"z": 1.1471657505990167, "p_value": 0.25131312369048175}

Go further

Add the dead-time correction of a high count rate (n = m / (1 - m.tau)).
Keep the results in a control chart and flag 2-sigma drifts with stat zscore.

executed on MentDB

Proton therapy: energy layers with the Bragg-Kleeman rule

To cover a tumour from 12 to 16 cm water-equivalent depth, the script converts each layer depth into a beam energy with the Bragg-Kleeman range-energy rule, sorts the layers, and gives the energy step and the range straggling at the distal edge.

Pipeline

01 · loopEach layer depthatom get→02 · energyE = (R/alpha)^(1/p)math power→03 · sortLayers in energy ordermath order_double→04 · meanMean energymath avg

Script

admin

#1. Proton therapy: tumour from 12 to 16 cm water-equivalent depth, layers every 0.8 cm, deepest first;
-> "[depths]" "[16, 15.2, 14.4, 13.6, 12.8, 12]";
#2. Bragg-Kleeman range R = alpha.E^p in water, alpha 0.0022 cm, p 1.77: E = R/alpha ^ 1/p;
-> "[list]" "";
for (-> "[i]" 1) (<= [i] 6) (++ "[i]") {
	-> "[r]" (atom get (string replace (string replace [depths] "[" "") "]" "") [i] ",");
	-> "[e]" (math round (math power (/ [r] 0.0022) (/ 1 1.77)) 1);
	-> "[list]" (concat [list] (if (> [i] 1) {", "} {""}) [e]);
};
-> "[es]" (math order_double (concat "[" [list] "]"));
#3. Range straggling about 1.2 % of the range at the distal edge, in mm;
concat "layer energies sorted = " [es] " MeV, mean = " (math round (math avg [es]) 1) " MeV, energy step at the distal edge = " (math round (- (math power (/ 16 0.0022) (/ 1 1.77)) (math power (/ 15.2 0.0022) (/ 1 1.77))) 2) " MeV, distal straggling = " (math round (* 0.012 160) 1) " mm"

Result

mentdb

layer energies sorted = [129.2,134.0,138.6,143.2,147.6,152.0] MeV, mean = 140.8 MeV, energy step at the distal edge = 4.34 MeV, distal straggling = 1.9 mm

Go further

Weight the layers to get a flat spread-out Bragg peak with matrix solve on a depth-dose matrix.
Compare with the PSTAR CSDA ranges with num interpolate.