-
Notifications
You must be signed in to change notification settings - Fork 16
SQ2.50 Procedural Framework: Part 7a (204 overcounts)
The first arithmetic in the "Part 7" code appears to relate to the fourth of the four fundamental controls. In SQUID 2.50, the calculations are carried out on a worksheet-by-worksheet basis, after SampleData has already been broken out into its sample-prefix sheets (if pbGrpAll = TRUE). So the below code is repeated for each sheet.
Frw = plaFirstDatRw(0) --"first data-row" (after column-header row)
Lrw = plaLastDatRw(0) --"last data-row"
If (pbExtractAgeGroups = TRUE) AND (piOverCtCorrType > 0)
--The dual condition here is curious: I see no a priori reason why pbExtractAgeGroups
--must be TRUE (although it is undoubtedly how SQUID 2.50 operates here). I think we
--should experiment with retaining only the second part of the If clause.
--Ludwig's verbatim comment for the overcount-correction equation is:
--OverctCorr 4/6 = (204cpsTot)-BkrdCps-OverctCPS)/(206cpsTot-BkrdCps)
--The parentheses are mismatched (first closing one should be deleted).
--"cpsTot" refers to "total counts at mass", without background subtraction.
For m = Frw To Lrw
GetSpotGroupingInfo Cells(m, 1)
--This subroutine updates the references to 204cps, BackgroundCps, 206cps, Nscans
--and the corresponding count_time_sec values, and will not be documented. Partly
--because it contains a bug that invalidates the calculation of the spot-by-spot
--uncertainty associated with the 204-overcount correction, but only when the
--original machine-data file was XML (PD files are not affected).
FinalTerm = "=1E-32"
If piOverCtCorrType = 1 --select appropriate BIWEIGHT MEAN value
FinalTerm = "= ( ["total 204 cts/sec"] - ["Bkrd cts/sec"] -
StandardData!Pb204OverCts7corr ) /
( ["total 206 cts/sec"] - ["Bkrd cts/sec"] )"
ElseIf piOverCtCorrType = 2
FinalTerm = "=( ["total 204 cts/sec"] - ["Bkrd cts/sec"] -
StandardData!Pb204OverCts8corr ) /
( ["total 206 cts/sec"] - ["Bkrd cts/sec"] )"
End If
PlaceFormulae FinalTerm, m, Pb46col --overwrites "raw" ["204/206"] with new formula
--As mentioned above, Squid3 should "archive" the original ["204/206"] under a
--different name (e.g. ["Raw 204/206"]), prior to this overwrite.
The next step is to calculate the uncertainty associated with overcount-corrected ["204/206"], but I believe Ludwig's arithmetic contains mathematical errors. I have sought to prove this, formulate an alternative, and demonstrate its correctness, at https://github.com/CIRDLES/Squid/issues/348. That page covers the issues in detail. The below code compiles both of the candidate expressions: terms defined by Ludwig's SQUID 2.50 math are suffixed 'KL'; my proposed revisions lack that suffix. The unclosed For loop continues:
If piNscans > 0 --trivial; I don't think the code can run this far otherwise
--Ludwig labels Term1KL as "var tot204cts":
Term1KL = ["total 204 cts/sec"] / ( piNscans * count_time_sec[204] )
--Reformulated correctly:
Term1 = ( ABS( ["total 204 cts/sec"] - ["Bkrd cts/sec"] ) /
( piNscans * count_time_sec[204] ) )
--Ludwig labels Term2KL as "var xs 204cps":
Term2KL = (StandardData!Pb204OverCts7corrEr)^2
--But hard-coding the value StandardData!Pb204OverCts7corrEr is incorrect:
--of course, the value used should depend on piOverCtCorrType. The following
--If clause is my own insertion, to address the bug.
If piOverCtCorrType = 1 --OverCts7corrEr
Term2 = (StandardData!Pb204OverCts7corrEr)^2
ElseIf piOverCtCorrType = 2 --OverCts8corrEr
Term2 = (StandardData!Pb204OverCts8corrEr)^2
End If
--Ludwig labels Term3KL as "var bkrdcps". Note that the below reference to
--["204/206"] targets the overcount-corrected value, not the "raw" value.
Term3KL = ( 1 - ["204/206"] )^2 * ["Bkrd cts/sec"] /
( piNscans * count_time_sec[Bkrd] )
--Reformulated correctly, Term3 depends on the value of piOverCtCorrType:
If piOverCtCorrType = 1 --OverCts7corr
Term3 = ( ( ["total 204 cts/sec"] - ["Bkrd cts/sec"] -
StandardData!Pb204OverCts7corr )^2 /
( piNscans * count_time_sec[206] *
ABS( ["total 206 cts/sec"] - ["Bkrd cts/sec"] ) )
ElseIf piOverCtCorrType = 2 --OverCts8corr
Term3 = ( ( ["total 204 cts/sec"] - ["Bkrd cts/sec"] -
StandardData!Pb204OverCts8corr )^2 /
( piNscans * count_time_sec[206] *
ABS( ["total 206 cts/sec"] - ["Bkrd cts/sec"] ) )
End If
--Ludwig does not label his Term4KL, which approximates Net206cps:
Term4KL = ( ["total 206 cts/sec"] - ["Bkrd cts/sec"] )
--Reformulated correctly, Term4 depends on the value of piOverCtCorrType:
If piOverCtCorrType = 1 --OverCts7corr
Term4 = ( ABS( ["total 204 cts/sec"] - ["Bkrd cts/sec"] -
StandardData!Pb204OverCts7corr ) )
ElseIf piOverCtCorrType = 2 --OverCts8corr
Term4 = ( ABS( ["total 204 cts/sec"] - ["Bkrd cts/sec"] -
StandardData!Pb204OverCts8corr ) )
End If
--Ludwig's Term5KL defines the expression for percentage uncertainty. Again,
--his reference to ["204/206"] targets the overcount-corrected (not "raw")
--value.
Term5KL = 100 * SQRT( Term1KL + Term2KL + Term3KL ) / Term4KL /
ABS( ["204/206"] )
--The reformulated version of Term5 makes no reference to the associated value;
--this is because the reformulated Terms1-4 are now fractional:
Term5 = 100 * SQRT( Term1 + Term2 + Term3 ) / Term4
Set Ce = Range(fsAddr(m, 1 + piPb46col)) --selects cell to the right of overcount-
--corrected ["204/206"], in order to place its %err value. As mentioned above,
--Squid3 should "archive" the original ["204/206 %err"] under a different name (e.g.
--["Raw 204/206 %err"]), prior to this overwrite.
--Finally, place the chosen expression for overcount-corrected ["204/206 %err"]:
--Ce.Formula = "=" & fsFcell(Term5KL, m) --LUDWIG expression, commented out
Ce.Formula = "=" & fsFcell(Term5, m) --REVISED expression, used in Squid3
--The next If tries to impose a range of 1e-9 < %err < 9999, but the variable
--"temp" appears unused in this context, so I am not sure it serves any purpose.
If IsNumeric(Ce)
temp = fvMax(0.000000001, fvMin(Ce, 9999))
End If
End If --piNscans > 0
Next m
Application.Calculate --refreshes all the arithmetic, now that the original ["204/206"]
--and its %err have been supplanted by their overcount-corrected analogues
--Finally, Ludwig amends the heading of the ["204/206"] column by adding the prefix
--"overct corr.", but I do NOT recommend this for Squid3. Note also that SQUID 2.50
--does not retain the measured ["204/206"] and its %err; the expressions for the
--overcount-corrected value and %err are written straight over the top. I really don't
--like that; as outlined above, I would much prefer instead that the "raw/measured"
--["204/206"] and its %err were archived unused to facilitate pre/post-correction
--comparisons, and that the overcount-corrected value and %err assume the "true" place
--of ["204/206"] and its %err in all the relevant arithmetic.
--But Ludwig does add a "hover-note" to the header, which I do like. It can be
--para-phrased as:
With Cells(Frw - 1, piPb46col)
If piOverCtCorrType = 1 --OverCts7corr
Note .Row, .Column, "Corrected for excess 204 counts inferred
from average 207Pb-corrected Standard spots"
ElseIf piOverCtCorrType = 2 --OverCts8corr
Note .Row, .Column, "Corrected for excess 204 counts inferred
from average 208Pb-corrected Standard spots"
End If
End With
End If --(pbExtractAgeGroups = TRUE) AND (piOverCtCorrType > 0), although see opening
--of If clause for the probable irrelevance of (pbExtractAgeGroups = TRUE)
As far as I can tell, this concludes the row-by-row implementation of the 204-overcount correction of ["204/206"] (when invoked) on Grouped-sample sheets. This first section of the procedure also introduces the concept of 'column-swapping' at isotopic-ration-column level, whereby a measured ["NUM/DEN"] and its %err can be substituted by user-defined arithmetical expression. This is relevant because SQUID 2.50 permits user-defined column-swapping as a Task Editor (i.e. Expression Manager in Squid3) feature. So the operation of "overcount-correcting" the ratio ["204/206"] can be considered an application of this feature that is sufficiently 'mainstream' that the U-Th-Pb data-reduction routine handles it (when invoked) as a matter of course. The extent to which Squid3 should support broader functionality (i.e. on any 'ratio of interest', according to user-defined arithmetic) remains an open question.