Index: trunk/tools/eam/teleff/cammask.pro
===================================================================
--- trunk/tools/eam/teleff/cammask.pro	(revision 41719)
+++ trunk/tools/eam/teleff/cammask.pro	(revision 41720)
@@ -92,4 +92,7 @@
  sprintf output "%s/%s.%s.cam.stf" $OUTDIR $expname $OUTVERSION
  write -fits DATA $output $OUT_FIELDS
+
+ # save this name for warp and diff analysis, if desired
+ $CAMSOURCE = $output
 end
 
Index: trunk/tools/eam/teleff/camstats.sh
===================================================================
--- trunk/tools/eam/teleff/camstats.sh	(revision 41719)
+++ trunk/tools/eam/teleff/camstats.sh	(revision 41720)
@@ -1,14 +1,19 @@
 
-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,97 +25,163 @@
    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,35 +192,42 @@
 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,5 +312,5 @@
 end
 
-macro load.camstats.glob
+macro load.camstats.glob.old
  if ($0 != 2)
    echo "USAGE: load.camstats.glob (glob)"
@@ -265,2 +343,25 @@
 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
Index: trunk/tools/eam/teleff/diffmask.pro
===================================================================
--- trunk/tools/eam/teleff/diffmask.pro	(revision 41719)
+++ trunk/tools/eam/teleff/diffmask.pro	(revision 41720)
@@ -13,4 +13,6 @@
 if (not($?MATCH_RADIUS)) set MATCH_RADIUS = 2.0
 
+if (not($?USE_CAMSOURCE)) set USE_CAMSOURCE = 0
+
 $SMF_FIELDS_DIF = X_PSF Y_PSF RA_PSF DEC_PSF CAL_PSF_MAG PSF_QF PSF_QF_PERFECT
 $REF_FIELDS_DIF = Rref Dref Mref 
@@ -50,5 +52,14 @@
   # load the reference source table
   # how to choose the name of the filter?
-  data $refsrcs; read -fits DATA $REF_FIELDS_DIF
+  data $refsrcs; 
+  if ($USE_CAMSOURCE) 
+   read -fits DATA Rref_out Dref_out Mref_out
+   set Rref = Rref_out
+   set Dref = Dref_out
+   set Mref = Mref_out
+  else
+    read -fits DATA $REF_FIELDS_DIF
+  end
+
   set Nref = ramp(Rref)
   if ($VERBOSE) echo "loaded Rref[] sources from reference $refsrcs"
Index: trunk/tools/eam/teleff/teleff.pro
===================================================================
--- trunk/tools/eam/teleff/teleff.pro	(revision 41719)
+++ trunk/tools/eam/teleff/teleff.pro	(revision 41720)
@@ -70,12 +70,17 @@
 
   # use a larger match radius for the asteroids
-  $MATCH_RADIUS = 3.0
+  # MATCH_RADIUS is in pixels : this should probably be even larger
+  $MATCH_RADIUS = 8.0
   $OUTVERSION = MPC
   $COPYMASKS = 0
 
-  make.truth.sources gpc1 $expname $mpcfile MPC
+  make.truth.sources gpc1 $expname $mpcfile MPC; # sets the value of REFSOURCE 
   get.camera.eff     gpc1 $expname EU.nightlyscience $REFSOURCE
   get.warp.eff       gpc1 $expname $camID            $REFSOURCE
   get.diff.eff       gpc1 $expname $warpID           $REFSOURCE
+  # load data from cam stage output:
+  $USE_CAMSOURCE = 1
+  get.warp.eff       gpc1 $expname $camID            $CAMSOURCE
+  get.diff.eff       gpc1 $expname $warpID           $CAMSOURCE
 
   return
Index: trunk/tools/eam/teleff/warpmask.pro
===================================================================
--- trunk/tools/eam/teleff/warpmask.pro	(revision 41719)
+++ trunk/tools/eam/teleff/warpmask.pro	(revision 41720)
@@ -5,4 +5,6 @@
 $OUTDIR     = .
 $OUTVERSION = XXX
+
+if (not($?USE_CAMSOURCE)) set USE_CAMSOURCE = 0
 
 $REF_FIELDS_WRP = Rref Dref Mref 
@@ -40,7 +42,15 @@
   # load the reference source table
   # how to choose the name of the filter?
-  data $refsrcs
   # NOTE: Mref is chosen to match filter for this exposure
-  data $refsrcs; read -fits DATA $REF_FIELDS_WRP
+  data $refsrcs                 ;
+  if ($USE_CAMSOURCE) 
+   read -fits DATA Rref_out Dref_out Mref_out
+   set Rref = Rref_out
+   set Dref = Dref_out
+   set Mref = Mref_out
+  else
+    read -fits DATA $REF_FIELDS_WRP
+  end
+
   set Nref = ramp(Rref)
   if ($VERBOSE) echo "loaded Rref[] sources from reference $refsrcs"
Index: trunk/tools/eam/teleff/warpstats.sh
===================================================================
--- trunk/tools/eam/teleff/warpstats.sh	(revision 41719)
+++ trunk/tools/eam/teleff/warpstats.sh	(revision 41720)
@@ -1,14 +1,19 @@
 
-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,80 +25,60 @@
    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,19 +86,19 @@
  $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
-