Skip to content

SQ2.50 Procedural Framework: Part 7a (204 overcounts)

sbodorkos edited this page Jul 22, 2019 · 2 revisions

Procedure

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.

Clone this wiki locally