Changeset 41720 for trunk/tools/eam/teleff/warpstats.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/warpstats.sh (modified) (3 diffs)

trunk/tools/eam/teleff/warpstats.sh

@@ -1 +1 @@
 
-macro warpstats.glob
+macro load.warpstats
  if ($0 != 2)
-   echo "USAGE: warpstats (glob)"
+   echo "USAGE: load.warpstats (glob)"
    break
  end
  
- delete -q Rref Dref Xref Yref Q1ref
+ # accumulate these fields for comparison plots
+ $FIELDS = 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.warpstats
    if ($i == 0)
-     set magFound_sum = magFound
      set numTotal_sum = numTotal
      set numPoss_sum  = numPoss
    else
+     set numTotal_sum = numTotal + numTotal_sum
      set numPoss_sum  = numPoss  + numPoss_sum
-     set numTotal_sum = numTotal + numTotal_sum
    end
-   concat RrefSS Rref
-   concat DrefSS Dref
-   concat XrefSS Xref
-   concat YrefSS Yref
-   concat Q1refSS Q1ref
+   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 fracPoss  = numPoss_sum  / numTotal_sum
+ set numPoss  = numPoss_sum
+ delete numTotal_sum numPoss_sum
 
- plot.warpstats
-
- # plot.warpsky
+ set fracPoss  = (numTotal > 0) ? numPoss  / numTotal : zero(numTotal)
 end
 
-macro warpstats.single
- if ($0 != 2)
-   echo "USAGE: warpstats (file)"
-   break
- end
- 
- data $1
- read -fits DATA -list-fields -q
- list tfields -join allfields
- read -fits DATA $allfields
+# assumes we have loaded values from the stf tables
+# magBin, numTotal, numFound, numPoss, numGood
+macro plot.warpstats.fracs
 
- # assumes data has been loaded 
- # results in magFound, fracFound, numFound
- calc.warpstats
-
- plot.warpstats
-end
-
-macro plot.warpstats
  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
+
+ # 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 5.5 +labelpady 5.0
- lim magFound -0.02 1.02; box $BSTY -ticks 1110;
- plot -x hist magFound fracPoss  -lw 3 -c blue70
- label -x mag -y "detectable fraction"
- subset tmp = fracPoss where (magFound > 15) && (magFound < 20); vstat -q tmp; $FracBrightPoss = $MEDIAN
+ lim magBin -0.02 1.02; box $BSTY -ticks 1110;
+
+ plot -x hist magBin fracPoss  -lw 3 -c blue70
+
+ subset tmp = fracPoss  where (magBin > 15) && (magBin < 20); vstat -q tmp; $FracBrightPoss  = $MEDIAN
  line -c blue -lt dot -lw 2 100 $FracBrightPoss to 0 $FracBrightPoss
- fprintf "%.1f-pct of sources are unmasked" {100*$FracBrightPoss}
 
+ fprintf "%.1f-pct of bright sources were detectable (unmasked)" {100*$FracBrightPoss}
+
+ label -x mag -y "detected fraction"
+
+ # plot the total number of detected sources, using the right axis
  section b 0 0 1 1
- lim magFound numTotal; box $BSTY -ticks 1011 -labels 0001; plot -x hist magFound numTotal -c grey70 -lw 2
- label +y "reference catalog numbers"
+ lim magBin numTotal; box $BSTY -ticks 1011 -labels 0001
+
+ plot -x hist magBin numTotal -c grey80 -lw 2
+ label +y "detected number"
+
+ $Xlegend = 0.02; $dXlegend = 0.03
+ $Ylegend = 0.95; $dYlegend = 0.02
+ line -frac $Xlegend $Ylegend to {$Xlegend + $dXlegend} $Ylegend -c blue  -lw 3; textline -frac -justify 5 {$Xlegend + $dXlegend + 0.01} $Ylegend "unmasked source positions (psfqf > 0.85)"; $Ylegend -= $dYlegend
 end
 
-macro plot.warpsky
-
-  dev -n sky
-  vstat -q Rref; $Ro = $MEDIAN
-  vstat -q Dref; $Do = $MEDIAN; $PlotRadius = 0.6*($MAX - $MIN)
-  region -n sky $Ro $Do $PlotRadius
-  cplot -pt box -sz 0.3 -op 0.1 -c black Rref Dref
-
-  # create t 0 360 0.01
-  # set dR = 3*dcos(t)/2.0
-  # set dD = 3*dsin(t)/2.0
-  # set Dc = $Do + dD
-  # set Rc = $Ro + dR/dcos(Dc)
-  # cplot -x line Rc Dc -c blue -lw 2
-end
-
-
+# relies on the following vectors: Mref_out, Q1ref_out, PSF_QF_out, X_PSF_out
 macro calc.warpstats
  $Mmin = 10
@@ -101 +86 @@
  $Mdel = 0.25
  
- delete -q magFound numTotal numPoss fracPoss
+ delete -q  magBin numTotal numPoss
  for mag $Mmin $Mmax $Mdel
-   concat $mag magFound
+   concat $mag magBin
 
-   # what fraction of the sources were actually detected?
-   set keepmag = (MrefSS > $mag - 0.5*$Mdel) && (MrefSS <= $mag + 0.5*$Mdel)
-   subset XinBin = XrefSS if keepmag
-   concat XinBin[] numTotal
+   # sources in mag bin
+   set keepSrc = (Mref_out > $mag - 0.5*$Mdel) && (Mref_out <= $mag + 0.5*$Mdel)
+   subset tmp = Mref_out if keepSrc
+   concat tmp[] numTotal
 
    # what fraction of the sources could have been detected (mask fraction > 0.85 at position)?
-   set keepgood = keepmag && (Q1refSS > 0.85)
-   subset XinBinGood = XrefSS if keepgood
-   concat {XinBinGood[]/XinBin[]} fracPoss
-   concat XinBinGood[] numPoss
+   set keepMsk = keepSrc && (Q1ref_out > 0.85)
+   subset tmp = Mref_out if keepMsk
+   concat tmp[] numPoss
  end
+
+ delete -q tmp keepSrc keepMsk keepDet
 end
-