Changeset 41720 for trunk/tools/eam/teleff/camstats.sh

Timestamp:

Jul 14, 2021, 3:30:26 PM (5 years ago)

Author:

eugene

Message:

improved plots, option to restrict to camera-confirmed sources

File:

• trunk/tools/eam/teleff/camstats.sh (modified) (5 diffs)

trunk/tools/eam/teleff/camstats.sh

@@ -1 +1 @@
 
-macro camstats.glob
+macro load.camstats
  if ($0 != 2)
-   echo "USAGE: camstats (glob)"
+   echo "USAGE: load.camstats (glob)"
    break
  end
  
- delete -q Rref Dref Xref Yref Q1ref Xpsf Ypsf Q1psf
+ # accumulate these fields for comparison plots
+ $FIELDS = RA_PSF DEC_PSF CAL_PSF_MAG PSF_QF Rref Dref Mref
+
+ # delete accumulation vectors
+ delete -q $FIELDS
 
  list name -glob $1
  for i 0 $name:n
 
+   echo "reading from $name:$i"
    data $name:$i
    read -fits DATA -sizes
@@ -20 +25 @@
    calc.camstats
    if ($i == 0)
-     set magFound_sum = magFound
+     set numTotal_sum = numTotal
      set numFound_sum = numFound
-     set numTotal_sum = numTotal
      set numPoss_sum  = numPoss
-     set numOK_sum    = numOK
+     set numGood_sum  = numGood
    else
-     set numOK_sum    = numOK    + numOK_sum
+     set numTotal_sum = numTotal + numTotal_sum
+     set numFound_sum = numFound + numFound_sum
      set numPoss_sum  = numPoss  + numPoss_sum
-     set numFound_sum = numFound + numFound_sum
-     set numTotal_sum = numTotal + numTotal_sum
+     set numGood_sum  = numGood  + numGood_sum
    end
-   concat RrefSS Rref
-   concat DrefSS Dref
-   concat XrefSS Xref
-   concat YrefSS Yref
-   concat Q1refSS Q1ref
-   concat X_PSF_found Xpsf
-   concat Y_PSF_found Ypsf
-   concat Q1_PSF_found Q1psf
+   foreach field $FIELDS
+     concat $field\_out $field
+   end
+
+   # cleanup by deleting the loaded vectors
+   delete -q $allfields
  end
    
- set magFound = magFound_sum
+ set numTotal = numTotal_sum
  set numFound = numFound_sum
- set numTotal = numTotal_sum
- set numOK    = numOK_sum
-
- set fracFound = numFound_sum / numTotal_sum
- set fracPoss  = numPoss_sum  / numTotal_sum
- set fracOK    = numOK_sum    / numTotal_sum
- set dfracFound = sqrt(numFound_sum) / numTotal_sum
- set dfracOK    = sqrt(numOK_sum)    / numTotal_sum
-
- break -auto off
- threshold -q -r magFound fracFound 0.25 -range 25 {magFound[]-1}
- break -auto on
- $MagLim = $threshold
-
- plot.camstats
-end
-
-macro camstats.single
- if ($0 != 2)
-   echo "USAGE: camstats (file)"
-   break
- end
- 
- data $1
- read -fits DATA -list-fields -q
- list tfields -join allfields
- read -fits DATA $allfields
-
- # assumes data has been loaded 
- # results in magFound, fracFound, numFound
- calc.camstats
-
- break -auto off
- threshold -q -r magFound fracFound 0.25  -range 25 {magFound[]-1}
- break -auto on
- $MagLim = $threshold
-
- plot.camstats
-end
-
-macro plot.camstats
+ set numPoss  = numPoss_sum
+ set numGood  = numGood_sum
+ delete numTotal_sum numFound_sum numPoss_sum numGood_sum
+
+ set fracFound = (numTotal > 0) ? numFound / numTotal : zero(numTotal)
+ set fracPoss  = (numTotal > 0) ? numPoss  / numTotal : zero(numTotal)
+ set fracGood  = (numTotal > 0) ? numGood  / numTotal : zero(numTotal)
+end
+
+if (not($?HISTBIN)) set HISTBIN = 0.01
+macro plot.camstats.resid
+
+  clear -s
+  resize 1800 1800
+  label -fn helvetica 18
+  $BST1 = -tickpad 0.2 -xpad 3.5 -labelpadx 3.0 -ypad 4.5 -labelpady 4.0 +ypad 0.5 +labelpady 0.0
+  $BST2 = -tickpad 0.2 -xpad 3.5 -labelpadx 3.0 -ypad 0.5 -labelpady 0.0 +ypad 4.5 +labelpady 4.0 -labels 1001
+  
+  section aM 0.00 0.66 0.65 0.33
+  section aR 0.00 0.33 0.65 0.33
+  section aD 0.00 0.00 0.65 0.33
+
+  section bM 0.65 0.66 0.35 0.33
+  section bR 0.65 0.33 0.35 0.33
+  section bD 0.65 0.00 0.35 0.33
+
+  set dR = 3600*(RA_PSF - Rref)*dcos(Dref)
+  set dD = 3600*(DEC_PSF - Dref)
+  set dM = CAL_PSF_MAG - Mref
+
+  set keepGood = (PSF_QF > 0.85)
+  set keepBright = keepGood && (Mref > 15) && (Mref < 17)
+
+  subset dRs = dR where keepBright
+  subset dDs = dD where keepBright
+  subset dMs = dM where keepBright
+
+  # allow mag and pos histograms to have different bin sizes?
+  histogram dRs NdRs -0.5 0.5 $HISTBIN -range dx
+  histogram dDs NdDs -0.5 0.5 $HISTBIN -range dx
+  histogram dMs NdMs -0.5 0.5 $HISTBIN -range dx
+
+  NdRs[0] = 0; NdRs[-1] = 0
+  NdDs[0] = 0; NdDs[-1] = 0
+  NdMs[0] = 0; NdMs[-1] = 0
+
+  $C3 = 0.0
+
+  peak -q dx NdRs; $C0 = $peakpos; $C2 = $peaknum; $C1 = 3*$HISTBIN; vgauss dx NdRs con NdRf
+  peak -q dx NdDs; $C0 = $peakpos; $C2 = $peaknum; $C1 = 3*$HISTBIN; vgauss dx NdDs con NdDf
+  peak -q dx NdMs; $C0 = $peakpos; $C2 = $peaknum; $C1 = 3*$HISTBIN; vgauss dx NdMs con NdMf
+
+  section aM; lim 13 23 -0.5 0.5; box $BST1
+  plot CAL_PSF_MAG dM -pt box -sz 0.5 -op 0.2 where keepGood
+  label -x mag -y "&sd&h mag"
+  
+  section aR; lim 13 23 -0.5 0.5; box $BST1
+  plot CAL_PSF_MAG dR -pt box -sz 0.5 -op 0.2 where keepGood
+  label -x mag -y "&sd&h R.A."
+  
+  section aD; lim 13 23 -0.5 0.5; box $BST1
+  plot CAL_PSF_MAG dD -pt box -sz 0.5 -op 0.2 where keepGood
+  label -x mag -y "&sd&h Dec"
+  
+  section bM; lim dx NdMs; box $BST2
+  plot dx NdMs -x hist -lw 3
+  plot dx NdMf -x line -c red -lw 2
+  label -x "&sd&h mag" +y "N (15 < mag < 17)"
+  
+  section bR; lim dx NdRs; box $BST2
+  plot dx NdRs -x hist -lw 3
+  plot dx NdRf -x line -c red -lw 2
+  label -x "&sd&h R.A" +y "N (15 < mag < 17)"
+  
+  section bD; lim dx NdDs; box $BST2
+  plot dx NdDs -x hist -lw 3
+  plot dx NdDf -x line -c red -lw 2
+  label -x "&sd&h Dec" +y "N (15 < mag < 17)"
+end
+
+# assumes we have loaded values from the stf tables
+# magBin, numTotal, numFound, numPoss, numGood
+macro plot.camstats.fracs
+
  clear -s
  resize 1800 1200
  label -fn helvetica 18
+ $BSTY = -tickpad 0.2 -xpad 3.5 -labelpadx 3.0 -ypad 4.5 -labelpady 4.0 +ypad 4.5 +labelpady 4.0
+
+ ## calculate the rough detection threshold (50%)
+ # find the bin with the most detections
+ peak -q magBin numFound
+
+ # search backwards along the fracFound curve until it crosses 50%
+ break -auto off
+ threshold -q -r magBin fracFound 0.5 -range $peaknum {magBin[]-1}
+ break -auto on
+ $MagLim = $threshold
+
+ # fraction error bars
+ set dfracFound = sqrt(numFound) / numTotal
+ set dfracGood  = sqrt(numGood)  / numTotal
+
+ # plot the detection fractions, using left axis
  section a 0 0 1 1
- $BSTY = -tickpad 0.2 -xpad 3.5 -labelpadx 3.0 -ypad 4.5 -labelpady 4.0 +ypad 4.5 +labelpady 4.0
- lim magFound -0.02 1.02; box $BSTY -ticks 1110;
- plot -x hist magFound fracFound -lw 3 -c red70
- plot -x hist magFound fracPoss  -lw 3 -c blue70
- plot -x hist magFound fracOK    -lw 3 -c red
+ lim magBin -0.02 1.02; box $BSTY -ticks 1110;
+
+ plot -x hist magBin fracFound -lw 3 -c red70
+ plot -x hist magBin fracPoss  -lw 3 -c blue70
+ plot -x hist magBin fracGood  -lw 3 -c red
+
+ subset tmp = fracGood  where (magBin > 15) && (magBin < 20); vstat -q tmp; $FracBrightGood  = $MEDIAN
+ subset tmp = fracPoss  where (magBin > 15) && (magBin < 20); vstat -q tmp; $FracBrightPoss  = $MEDIAN
+ subset tmp = fracFound where (magBin > 15) && (magBin < 20); vstat -q tmp; $FracBrightFound = $MEDIAN
+
+ line -c red  -lt dot -lw 2 100 $FracBrightGood to 0 $FracBrightGood
+ line -c blue -lt dot -lw 2 100 $FracBrightPoss to 0 $FracBrightPoss
+
+ fprintf "%.1f-pct of all bright sources were detected unmasked" {100*$FracBrightGood} 
+ fprintf "%.1f-pct of unmasked bright sources were found" {100*$FracBrightFound / $FracBrightPoss}
+ fprintf "%.1f-pct of unmasked bright sources were detected unmasked" {100*$FracBrightGood / $FracBrightPoss}
+ fprintf "mag limit %.1f" $MagLim
+
  label -x mag -y "detected fraction"
- subset tmp = fracOK    where (magFound > 15) && (magFound < 20); vstat -q tmp; $FracBrightOK    = $MEDIAN
- subset tmp = fracPoss  where (magFound > 15) && (magFound < 20); vstat -q tmp; $FracBrightPoss  = $MEDIAN
- subset tmp = fracFound where (magFound > 15) && (magFound < 20); vstat -q tmp; $FracBrightFound = $MEDIAN
- line -c red  -lt dot -lw 2 100 $FracBrightOK to 0 $FracBrightOK
- line -c blue -lt dot -lw 2 100 $FracBrightPoss to 0 $FracBrightPoss
- fprintf "%.1f-pct of all bright sources were detected unmasked" {100*$FracBrightOK} 
- fprintf "%.1f-pct of unmasked bright sources were found" {100*$FracBrightFound / $FracBrightPoss}
- fprintf "%.1f-pct of unmasked bright sources were detected unmasked" {100*$FracBrightOK / $FracBrightPoss}
- fprintf "mag limit %.1f" $MagLim
-
+
+ # plot the total number of detected sources, using the right axis
  section b 0 0 1 1
- lim magFound numFound; box $BSTY -ticks 1011 -labels 0001; plot -x hist magFound numFound -c grey70 -lw 2
+ lim magBin numTotal; box $BSTY -ticks 1011 -labels 0001
+
+ plot -x hist magBin numFound -c grey50 -lw 2
+ plot -x hist magBin numTotal -c grey80 -lw 2
+
  line -c red -lt dot -lw 2 $MagLim -100 to $MagLim 10000
  label +y "detected number"
 
  vstat -q numFound; $Nfound = $TOTAL
- vstat -q numOK;    $Nok    = $TOTAL
- 
- fprintf "%d total sources detected, %d unmasked" $Nfound $Nok
+ vstat -q numGood;  $Ngood  = $TOTAL
+ 
+ fprintf "%d total sources detected, %d unmasked" $Nfound $Ngood
 
  $Xlegend = 0.02; $dXlegend = 0.03
@@ -121 +192 @@
 end
 
+# generate histograms as a function of magnitude for the following:
+# number of sources (numTotal)
+# number of sources which were not masked, possibly detected (numPoss)
+# number detected (numFound)
+# number of detections with good quality (numGood)
+
+# relies on the following vectors: Mref_out, Q1ref_out, PSF_QF_out, X_PSF_out
 macro calc.camstats
  $Mmin = 10
  $Mmax = 25
- $Mdel = 0.5
- 
- delete -q fracFound magFound numFound numTotal numPoss fracPoss numOK fracOK
+ $Mdel = 0.25
+ 
+ delete -q  magBin numTotal numPoss numFound numGood
  for mag $Mmin $Mmax $Mdel
-   concat $mag magFound
-
-   # what fraction of the sources were actually detected?
-   set keep = (MrefSS > $mag - 0.5*$Mdel) && (MrefSS <= $mag + 0.5*$Mdel)
-   # X_PSF_found is NAN if the source is not found (detected in smf), finite otherwise
-   subset XinBin = X_PSF_found if keep
-   vstat -q XinBin
-   concat {$NPTS/XinBin[]} fracFound
-   concat XinBin[] numTotal
+   concat $mag magBin
+
+   # sources in mag bin
+   set keepSrc = (Mref_out > $mag - 0.5*$Mdel) && (Mref_out <= $mag + 0.5*$Mdel)
+   subset tmp = X_PSF_out if keepSrc
+   concat tmp[] numTotal
+
+   # NOTE: X_PSF_out is NAN if the source is not found (detected in smf), finite otherwise
+   # vstat ignores NAN values, so NPTS is count of non-NAN values
+   vstat -q tmp
    concat $NPTS numFound
 
+   # what fraction of the sources could have been detected (mask fraction > 0.85 at position)?
+   set keepMsk = keepSrc && (Q1ref_out > 0.85)
+   subset tmp = X_PSF_out if keepMsk
+   concat tmp[] numPoss
+
    # what fraction of the sources were actually detected with good quality?
-   set keep = (MrefSS > $mag - 0.5*$Mdel) && (MrefSS <= $mag + 0.5*$Mdel) && (Q1_PSF_found > 0.85) && not(isnan(Q1_PSF_found))
-   # X_PSF_found is NAN if the source is not found (detected in smf), finite otherwise
-   subset XinBinOK = X_PSF_found if keep
-   concat {XinBinOK[]/XinBin[]} fracOK
-   concat XinBinOK[] numOK
-
-   # what fraction of the sources could have been detected (mask fraction > 0.85 at position)?
-   set keep = (MrefSS > $mag - 0.5*$Mdel) && (MrefSS <= $mag + 0.5*$Mdel) && (Q1refSS > 0.85)
-   subset XinBinGood = X_PSF_found if keep
-   concat {XinBinGood[]/XinBin[]} fracPoss
-   concat XinBinGood[] numPoss
- end
+   set keepDet = keepSrc && (PSF_QF_out > 0.85) && not(isnan(PSF_QF_out))
+   subset tmp = X_PSF_out if keepDet
+   concat tmp[] numGood
+ end
+
+ delete -q tmp keepSrc keepMsk keepDet
 end
 
@@ -234 +312 @@
 end
 
-macro load.camstats.glob
+macro load.camstats.glob.old
  if ($0 != 2)
    echo "USAGE: load.camstats.glob (glob)"
@@ -265 +343 @@
 end
 
+# relies on the following vectors: Mref_out, Q1ref_out, PSF_QF_out, X_PSF_out
+macro test.camstats
+ $Mmin = 10
+ $Mmax = 25
+ $Mdel = 0.25
+ 
+ delete -q  magBin numTotal numFound
+ for mag $Mmin $Mmax $Mdel
+   concat $mag magBin
+
+   # sources in mag bin
+   set keepSrc = (MrefDiff_onchp > $mag - 0.5*$Mdel) && (MrefDiff_onchp <= $mag + 0.5*$Mdel)
+   subset tmp = RpsfDiff_onchp if keepSrc
+   concat tmp[] numTotal
+
+   # NOTE: X_PSF_out is NAN if the source is not found (detected in smf), finite otherwise
+   # vstat ignores NAN values, so NPTS is count of non-NAN values
+   vstat -q tmp
+   concat $NPTS numFound
+ end
+
+ delete -q tmp keepSrc
+end