Index: /branches/eam_branches/ipp-20211108/ippTests/ippCheck/.dvo
===================================================================
--- /branches/eam_branches/ipp-20211108/ippTests/ippCheck/.dvo	(revision 41941)
+++ /branches/eam_branches/ipp-20211108/ippTests/ippCheck/.dvo	(revision 41941)
@@ -0,0 +1,1 @@
+echo $MODULES:0
Index: /branches/eam_branches/ipp-20211108/ippTests/ippCheck/.mana
===================================================================
--- /branches/eam_branches/ipp-20211108/ippTests/ippCheck/.mana	(revision 41941)
+++ /branches/eam_branches/ipp-20211108/ippTests/ippCheck/.mana	(revision 41941)
@@ -0,0 +1,12 @@
+echo $MODULES:n
+echo $MODULES:
+echo $MODULES:0
+!ls $MODULES:0
+!ls $MODULES:0/plots
+!less $MODULES:0
+!less $MODULES:0/plo
+!less $MODULES:0/plots
+load plots
+which load
+module
+input
Index: /branches/eam_branches/ipp-20211108/ippTests/ippCheck/.pantasks
===================================================================
--- /branches/eam_branches/ipp-20211108/ippTests/ippCheck/.pantasks	(revision 41941)
+++ /branches/eam_branches/ipp-20211108/ippTests/ippCheck/.pantasks	(revision 41941)
@@ -0,0 +1,1 @@
+echo $MODULES:0
Index: /branches/eam_branches/ipp-20211108/ippTests/ippCheck/Makefile.am
===================================================================
--- /branches/eam_branches/ipp-20211108/ippTests/ippCheck/Makefile.am	(revision 41941)
+++ /branches/eam_branches/ipp-20211108/ippTests/ippCheck/Makefile.am	(revision 41941)
@@ -0,0 +1,18 @@
+module_files = \
+  cammask.pro   /
+  camstats.pro  /
+  diffmask.pro  /
+  truth.pro     /
+  warpmask.pro  /
+  warpstats.pro
+
+dvo_dir = $(datadir)/dvo
+module_dir = $(dvo_dir)/modules
+ippchecks_DATA = $(module_files)
+
+install-data-hook:
+	chmod 0755 $(dvo_dir) $(module_dir)
+
+EXTRA_DIST = $(module_files)
+
+ACLOCAL_AMFLAGS = -I m4
Index: /branches/eam_branches/ipp-20211108/ippTests/ippCheck/autogen.sh
===================================================================
--- /branches/eam_branches/ipp-20211108/ippTests/ippCheck/autogen.sh	(revision 41941)
+++ /branches/eam_branches/ipp-20211108/ippTests/ippCheck/autogen.sh	(revision 41941)
@@ -0,0 +1,103 @@
+#!/bin/sh
+# Run this to generate all the initial makefiles, etc.
+
+srcdir=`dirname $0`
+test -z "$srcdir" && srcdir=.
+
+ORIGDIR=`pwd`
+cd $srcdir
+
+PROJECT=ippCheck
+TEST_TYPE=-f
+# change this to be a unique filename in the top level dir
+FILE=autogen.sh
+
+DIE=0
+
+LIBTOOLIZE=libtoolize
+ACLOCAL="aclocal $ACLOCAL_FLAGS"
+AUTOHEADER=autoheader
+AUTOMAKE=automake
+AUTOCONF=autoconf
+
+#($LIBTOOLIZE --version) < /dev/null > /dev/null 2>&1 || {
+#        echo
+#        echo "You must have $LIBTOOLIZE installed to compile $PROJECT."
+#        echo "Download the appropriate package for your distribution,"
+#        echo "or get the source tarball at http://ftp.gnu.org/gnu/libtool/"
+#        DIE=1
+#}
+
+($ACLOCAL --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have $ACLOCAL installed to compile $PROJECT."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/"
+        DIE=1
+}
+
+#($AUTOHEADER --version) < /dev/null > /dev/null 2>&1 || {
+#        echo
+#        echo "You must have $AUTOHEADER installed to compile $PROJECT."
+#        echo "Download the appropriate package for your distribution,"
+#        echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/"
+#        DIE=1
+#}
+
+($AUTOMAKE --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have $AUTOMAKE installed to compile $PROJECT."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/"
+        DIE=1
+}
+
+($AUTOCONF --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have $AUTOCONF installed to compile $PROJECT."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/"
+        DIE=1
+}
+
+if test "$DIE" -eq 1; then
+        exit 1
+fi
+
+test $TEST_TYPE $FILE || {
+        echo "You must run this script in the top-level $PROJECT directory"
+        exit 1
+}
+
+if test -z "$*"; then
+        echo "I am going to run ./configure with no arguments - if you wish "
+        echo "to pass any to it, please specify them on the $0 command line."
+fi
+
+#$LIBTOOLIZE --copy --force || echo "$LIBTOOLIZE failed"
+$ACLOCAL || echo "$ACLOCAL failed"
+#$AUTOHEADER || echo "$AUTOHEADER failed"
+$AUTOMAKE --add-missing --force-missing --copy || echo "$AUTOMAKE failed"
+$AUTOCONF || echo "$AUTOCONF failed"
+
+cd $ORIGDIR
+
+run_configure=true
+for arg in $*; do
+    case $arg in
+        --no-configure)
+            run_configure=false
+            ;;
+        *)
+            ;;
+    esac
+done
+
+if $run_configure; then
+    $srcdir/configure --enable-maintainer-mode "$@"
+    echo
+    echo "Now type 'make' to compile $PROJECT."
+else
+    echo
+    echo "Now run 'configure' and 'make' to compile $PROJECT."
+fi
Index: /branches/eam_branches/ipp-20211108/ippTests/ippCheck/cammask.pro
===================================================================
--- /branches/eam_branches/ipp-20211108/ippTests/ippCheck/cammask.pro	(revision 41941)
+++ /branches/eam_branches/ipp-20211108/ippTests/ippCheck/cammask.pro	(revision 41941)
@@ -0,0 +1,380 @@
+#!/usr/bin/env dvo
+
+$VERBOSE    = 1	       
+$COPYMASKS  = 0
+$OUTDIR     = .
+$OUTVERSION = XXX
+
+if (not($?MATCH_RADIUS)) set MATCH_RADIUS = 2.0
+
+$SMF_FIELDS_CAM = X_PSF Y_PSF RA_PSF DEC_PSF CAL_PSF_MAG PSF_QF PSF_QF_PERFECT
+$REF_FIELDS_CAM = Rref Dref Mref 
+$XTR_FIELDS_CAM = Xref Yref Nref Q1ref
+
+macro get.camera.eff
+ if ($0 != 5)
+  echo "USAGE: get.camera.eff (camera) (expname) (camLabel) (refsrcs.fits)"
+  echo "EX:    get.camera.eff gpc1 o9322g0253o EU.nightlyscience filename.fits"
+  break
+ end
+
+ local i dbname expname camlabel output
+ $dbname   = $1
+ $expname  = $2
+ $camlabel = $3
+ $refsrcs  = $4
+
+ $OUTDIR = $expname.stats
+ mkdir $OUTDIR
+
+ # I have a table of sources which may be in the exposure
+ # I need to generate a matched table of all sources in the smf tables
+ # AND the mask fraction (roughly psf_qf) for each detection 
+
+ # get the name of the smf file:
+ get.camera.metadata $dbname $expname $camlabel
+ $smffile = $camPath.smf     
+ $smfreal = `neb-locate -p $smffile`
+ 
+ get.chip.metadata $dbname $chipID
+ get.chip.mask.files $camPath
+ 
+ for i 0 $camMask:n
+   echo $chipName:$i $camMask:$i
+ end
+
+ # load the reference source table
+ # NOTE: Mref is chosen to match filter for this exposure
+ data $refsrcs; 
+ read -fits DATA -sizes
+ if ($table:Ny == 0)
+  echo "ERROR: no reference sources found in $refsrcs"
+  break
+ end
+ 
+ read -fits DATA $REF_FIELDS_CAM
+ set Nref = ramp(Rref)
+
+ # read in the PHU header
+ rd PHU $smfreal
+
+ # create the chipName list from the smffile and calculate the index numbers.
+ # NOTE: if a chip is missing from the smf, this may be inconsistent with chipMask:n
+ # XXX for now, we fail is this happens
+ list chipsRaw -x "ftable -list $smfreal | grep -n hdr | prcol 1 | sed s/.hdr//"
+ for i 0 $chipsRaw:n
+  list word -splitbychar : $chipsRaw:$i
+  $chipNseq:$i = $word:0 - 3; # this number depends on the number of leading tables 
+  if ("$chipName:$i" != "$word:1")
+    echo "ERROR: missing chip from smffile: " $chipName:$i
+    break
+  end
+ end
+
+  # delete any existing output vectors
+  foreach field $SMF_FIELDS_CAM $REF_FIELDS_CAM $XTR_FIELDS_CAM
+    delete -q $field\_out
+  end
+
+ # find matching entries in smf file 
+ for i 0 $chipName:n
+   get.chip.detect.stats $chipName:$i $chipNseq:$i
+   get.chip.mask.stats $camMask:$i
+
+   foreach field $SMF_FIELDS_CAM
+     concat $field\_found $field\_out
+     delete $field\_found
+   end
+   foreach field $REF_FIELDS_CAM $XTR_FIELDS_CAM
+     concat $field\_onchp $field\_out
+     delete $field\_onchp
+   end
+ end
+
+ $OUT_FIELDS = 
+ foreach field $SMF_FIELDS_CAM $REF_FIELDS_CAM $XTR_FIELDS_CAM
+   $OUT_FIELDS = $OUT_FIELDS $field\_out
+ end
+
+ 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
+
+# We can read the data by EXTNAME or sequence number. Reading by
+# EXTNAME is more robust against format changes but slow
+$USE_EXTNAME = 0
+
+macro get.chip.detect.stats
+  if ($0 != 3) 
+   echo "USAGE: get.chip.detect.stats (chipname) (chipnseq)"
+   echo "  assumes existence of $smfreal vectors (Rref, Dref, Mref) and buffer (PHU)"
+   echo "  matches known sources to detections in the smf file extensions for each chip"
+   echo "  calculates detection fraction as a function of magnitude"
+   echo "  also astrometric and photometric scatter vs magnitude for good detections"
+   break
+  end
+
+  local chipname chipnseq field NAXIS1 NAXIS2
+  $chipname = $1
+  $chipnseq = $2
+
+  # read in the CHIP header
+
+  if ($USE_EXTNAME) 
+    sprintf extname "%s.hdr" $chipname
+    rd -n $extname hdr $smfreal
+  else
+    rd -x $chipnseq hdr $smfreal
+  end
+
+  keyword hdr NAXIS1 NAXIS1
+  keyword hdr NAXIS2 NAXIS2
+
+  coords -q hdr -mosaic PHU -p       0       0; init.radec.range
+  coords -q hdr -mosaic PHU -p $NAXIS1       0; set.radec.range
+  coords -q hdr -mosaic PHU -p $NAXIS1 $NAXIS2; set.radec.range
+  coords -q hdr -mosaic PHU -p       0 $NAXIS2; set.radec.range
+
+  # generate a subset which could be on this chip
+  # do we need to grow the range at all?
+  set keep = (Rref >= $minR) && (Rref <= $maxR) && (Dref >= $minD) && (Dref <= $maxD) 
+  foreach field $REF_FIELDS_CAM Nref
+    subset $field\_tmp = $field if keep
+  end
+
+  # convert the subset coordinates to pixels
+  set Xref_tmp = Rref_tmp
+  set Yref_tmp = Dref_tmp
+  coords -q hdr -mosaic PHU -c Xref_tmp Yref_tmp
+
+  # generate a subset of sources which are within the boundaries of the chip
+  set keep = (Xref_tmp >= 0) && (Xref_tmp <= $NAXIS1) && (Yref_tmp >= 0) && (Yref_tmp <= $NAXIS2)
+  foreach field $REF_FIELDS_CAM Nref Xref Yref
+    subset $field\_onchp = $field\_tmp if keep
+    delete $field\_tmp
+  end
+
+  # now read the smf detections
+  data $smfreal
+
+  if ($USE_EXTNAME) 
+    sprintf extname "%s.psf" $chipname
+    read -fits $extname $SMF_FIELDS_CAM
+  else
+    read -fits {$chipnseq+1} -extnum $SMF_FIELDS_CAM
+  end
+
+  # match radius currently needs to be larger for MPC objects
+  match2d Xref_onchp Yref_onchp X_PSF Y_PSF $MATCH_RADIUS -closest
+  foreach field $SMF_FIELDS_CAM
+    reindex $field\_found = $field using index1 -keep-unmatched
+  end
+end
+
+macro get.chip.mask.stats
+  if ($0 != 2) 
+   echo "USAGE: get.chip.mask.stats (cammask)"
+   echo "  assumes existence of source chip positions (Xref_onchp, Yref_onchp)"
+   break
+  end
+
+  local cammask camreal
+  $cammask = $1
+
+  $camreal = `neb-locate -p $cammask`
+  rd -x 0 mk $camreal 
+
+  # I have a mask file, now I need to calculate something equiv to psfQf
+  set mkPerfect = (mk == 0) ? zero(mk) + 1 : zero(mk)
+
+  # I want to set an image based on the bit values but I fear
+  # the mask has been converted to a float when read
+  forcedphot mkPerfect Xref_onchp Yref_onchp -output Q1ref_onchp
+end
+
+### utilities
+macro init.radec.range
+  $minR = $RA
+  $minD = $DEC
+  $maxR = $RA
+  $maxD = $DEC
+end
+
+macro set.radec.range
+  $minR = min($minR, $RA)
+  $minD = min($minD, $DEC)
+  $maxR = max($maxR, $RA)
+  $maxD = max($maxD, $DEC)
+end
+
+# get the file names of the chip mask files.
+# these are actually the masks from the camera stage which
+# has the dynamic masking applied
+macro get.chip.mask.files
+ if ($0 != 2)
+  echo "USAGE: get.chip.mask.files (camPath)"
+  echo " uses chipName:0 - chipName:n"
+  echo " sets camMask:0  - camMask:n"
+  break
+ end
+
+ local pathbase i j mychip Nch
+ $pathbase = $1
+
+ # get the full list from nebulous (faster than getting each one at a time)
+ list tmpName -x "neb-ls $pathbase.%.mk.fits"
+
+ # the list of masks might not match the list of chips
+ $j = 0
+ for i 0 $chipName:n
+   list word -splitbychar . $tmpName:$j
+   $Nch = $word:n - 3
+   $mychip = $word:$Nch
+   if ($mychip != $chipName:$i) 
+    $camMask:$i = NONE
+    continue
+   end
+   $camMask:$i = $tmpName:$j
+   if ($VERBOSE) echo "found chipMask $camMask:$i"
+   if ($COPYMASKS) exec neb-copy $camMask:$i $OUTDIR
+   $j ++
+ end
+ $camMask:n = $chipName:n
+end
+
+macro get.chip.metadata
+ if ($0 != 3)
+  echo "USAGE: get.chip.metadata (camera) (chipID)"
+  echo "EX:    get.chip.metadata gpc1 NNNN"
+  echo "  sets camPath, chipID"
+  break
+ end
+
+ local dbname chipID
+ $dbname  = $1
+ $chipID  = $2
+
+ $command = chiptool -dbname $dbname -processedimfile -chip_id $chipID
+ 
+ queueload  myqueue -x "$command"
+ queue2book myqueue mybook
+ 
+ book npages mybook -var npages
+ 
+ if ($npages < 1) 
+  echo "ERROR: missing exposure : $command"
+  book delete mybook
+  return
+ end
+ 
+ local i
+ for i 0 $npages 
+   sprint pagename page.%03d $i
+   book getword mybook $pagename class_id -var chipName:$i
+ end
+ $chipName:n = $npages
+
+ book delete mybook
+end
+
+macro get.camera.metadata
+ if ($0 != 4)
+  echo "USAGE: get.camera.metadata (camera) (expname) (camLabel)"
+  echo "EX:    get.camera.metadata gpc1 o9322g0253o EU.nightlyscience"
+  echo "  sets camPath, chipID"
+  break
+ end
+
+ local dbname expname camlabel
+ $dbname   = $1
+ $expname  = $2
+ $camlabel = $3
+
+# need to generically handle the question of which version / label
+# this will be a problem for old exposures without a useful label (cleanup)
+# instead of label, we could specify the cam_id (which we might know)
+
+ $command = camtool -dbname $dbname -processedexp -exp_name $expname -label $camlabel
+#$command = camtool -dbname $dbname -processedexp -exp_name $expname
+ 
+ queueload  myqueue -x "$command"
+ queue2book myqueue mybook
+ 
+ $camPath = NONE
+ book npages mybook -var npages
+ 
+ if ($npages < 1) 
+  echo "ERROR: missing exposure : $command"
+  book delete mybook
+  next
+ end
+ 
+ if ($npages > 1) 
+  echo "ERROR: too many exposures : $command"
+  book delete mybook
+  next
+ end
+ 
+ sprint pagename page.%03d 0
+ book getword mybook $pagename path_base -var camPath
+ book getword mybook $pagename chip_id   -var chipID
+ book getword mybook $pagename cam_id    -var camID
+ book delete mybook
+end
+
+
+#### this test macro generates the match table without measuring the mask fraction
+#### XXX the output formats are not consistent with the above version
+
+macro get.camera.matches
+ if ($0 != 5)
+  echo "USAGE: get.camera.matches (smfreal) (refsrcs.fits) (outroot) (outversion)"
+  break
+ end
+
+ local i refsrcs 
+ $smfreal = $1
+ $refsrcs = $2
+ $outroot = $3
+ $OUTVERSION = $4
+
+ dirname $outroot -var OUTDIR
+ mkdir $OUTDIR
+
+ # load the reference source table
+ # NOTE: Mref is chosen to match filter for this exposure
+ data $refsrcs; read -fits DATA $REF_FIELDS_CAM
+ set Nref = ramp(Rref)
+
+ # read in the PHU header
+ rd PHU $smfreal
+
+ # create the chipName list from the smffile and calculate the index numbers
+ list chipsRaw -x "ftable -list $smfreal | grep -n hdr | prcol 1 | sed s/.hdr//"
+ for i 0 $chipsRaw:n
+  list word -splitbychar : $chipsRaw:$i
+  $chipNseq:$i = $word:0 - 3
+  $chipName:$i = $word:1
+ end
+ $chipName:n = $chipsRaw:n
+ $chipNseq:n = $chipsRaw:n
+
+ # find matching entries in smf file 
+ for i 0 $chipName:n
+   echo $chipName:$i
+   get.chip.detect.stats $chipName:$i $chipNseq:$i
+   # if we do not have masks, we have not measured Q1refSS, just set to 1.0
+   set Q1ref_onchp = zero(Rref_onchp) + 1.0
+
+ end
+
+ $REF_FIELDS_OUT = $REF_FIELDS Xref Yref Q1ref
+
+ sprintf output "%s.%s.%s.cam.stf" $outroot $chipName:$i $OUTVERSION
+ write -fits DATA $output RrefSS DrefSS MrefSS XrefSS YrefSS Q1refSS X_PSF_found Y_PSF_found M_PSF_found RA_PSF_found DEC_PSF_found Q1_PSF_found Q2_PSF_found 
+
+end
+
Index: /branches/eam_branches/ipp-20211108/ippTests/ippCheck/camstats.pro
===================================================================
--- /branches/eam_branches/ipp-20211108/ippTests/ippCheck/camstats.pro	(revision 41941)
+++ /branches/eam_branches/ipp-20211108/ippTests/ippCheck/camstats.pro	(revision 41941)
@@ -0,0 +1,405 @@
+
+macro load.camstats
+ if ($0 != 2)
+   echo "USAGE: load.camstats (glob)"
+   break
+ end
+ 
+ # 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
+   if ($table:Ny == 0) continue
+   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
+   if ($i == 0)
+     set numTotal_sum = numTotal
+     set numFound_sum = numFound
+     set numPoss_sum  = numPoss
+     set numGood_sum  = numGood
+   else
+     set numTotal_sum = numTotal + numTotal_sum
+     set numFound_sum = numFound + numFound_sum
+     set numPoss_sum  = numPoss  + numPoss_sum
+     set numGood_sum  = numGood  + numGood_sum
+   end
+   foreach field $FIELDS
+     concat $field\_out $field
+   end
+
+   # cleanup by deleting the loaded vectors
+   delete -q $allfields
+ end
+   
+ set numTotal = numTotal_sum
+ set numFound = numFound_sum
+ 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
+if (not($?BRIGHTMIN)) set BRIGHTMIN = 15.0
+if (not($?BRIGHTMAX)) set BRIGHTMAX = 17.0
+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 > $BRIGHTMIN) && (Mref < $BRIGHTMAX)
+
+  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 {-50*$HISTBIN} {50*$HISTBIN} $HISTBIN -range dx
+  histogram dDs NdDs {-50*$HISTBIN} {50*$HISTBIN} $HISTBIN -range dx
+  histogram dMs NdMs {-50*$HISTBIN} {50*$HISTBIN} $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 = 2*$HISTBIN; vgauss dx NdRs con NdRf; $SigmaRA  = abs($C1)
+  peak -q dx NdDs; $C0 = $peakpos; $C2 = $peaknum; $C1 = 2*$HISTBIN; vgauss dx NdDs con NdDf; $SigmaDec = abs($C1)
+  peak -q dx NdMs; $C0 = $peakpos; $C2 = $peaknum; $C1 = 2*$HISTBIN; vgauss dx NdMs con NdMf; $SigmaMag = abs($C1)
+
+  section aM; lim 13 23 {-50*$HISTBIN} {50*$HISTBIN}; 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 {-50*$HISTBIN} {50*$HISTBIN}; 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 {-50*$HISTBIN} {50*$HISTBIN}; 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 ($BRIGHTMIN < mag < $BRIGHTMAX)"
+  sprintf line "&ss&h (mag) = %.0f mmag" {1000*$SigmaMag}; textline -frac 0.02 0.95 -justify 5 "$line"
+  
+  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 ($BRIGHTMIN < mag < $BRIGHTMAX)"
+  sprintf line "&ss&h (mag) = %.0f mas" {1000*$SigmaRA}; textline -frac 0.02 0.95 -justify 5 "$line"
+  
+  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 ($BRIGHTMIN < mag < $BRIGHTMAX)"
+  sprintf line "&ss&h (Dec) = %.0f mas" {1000*$SigmaDec}; textline -frac 0.02 0.95 -justify 5 "$line"
+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 = 0.0
+ 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
+ lim magBin -0.02 1.02; box $BSTY -ticks 1110;
+
+ plot -x hist magBin fracPoss  -lw 3 -c blue70
+ plot -x hist magBin fracFound -lw 3 -c red70
+ plot -x hist magBin fracGood  -lw 3 -c red
+
+ 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
+ subset tmp = fracGood  where (magBin > 15) && (magBin < 20); vstat -q tmp; $FracBrightGood  = $MEDIAN
+
+ line -c blue70 -lt dot -lw 2 100 $FracBrightPoss  to 0 $FracBrightPoss
+ line -c red70  -lt dot -lw 2 100 $FracBrightFound to 0 $FracBrightFound
+ line -c red    -lt dot -lw 2 100 $FracBrightGood  to 0 $FracBrightGood
+
+ 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
+ sprintf line "mag limit: %.1f" $MagLim
+ textline -frac 0.95 0.95 -justify 3 "$line"
+
+ label -x mag -y "detected fraction"
+
+ # plot the total number of detected sources, using the right axis
+ section b 0 0 1 1
+ 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 numTotal; $Ntotal = $TOTAL
+ vstat -q numFound; $Nfound = $TOTAL
+ vstat -q numGood;  $Ngood  = $TOTAL
+ 
+ fprintf "%d of %d total sources detected, %d unmasked" $Nfound $Ntotal $Ngood
+
+ legend.init 0.02 0.95 0.03 0.02
+ sprintf line "unmasked source positions (psfqf > 0.85) [%.1f pct]" {100*$FracBrightPoss};  legend.line blue  3 solid "$line"
+ sprintf line "all sources detected [%.1f pct]"                     {100*$FracBrightFound}; legend.line red70 3 solid "$line"
+ sprintf line "sources detected with psfqf > 0.85 [%.1f pct]"	    {100*$FracBrightGood};  legend.line red   3 solid "$line"
+end
+
+macro legend.init
+ if ($0 != 5)
+  echo "USAGE: legend.init (X) (Y) (dX) (dY)"
+  break
+ end
+ 
+ $legend:X = $1; $legend:dX = $3
+ $legend:Y = $2; $legend:dY = $4
+end
+
+macro legend.line
+ if ($0 != 5)
+  echo "USAGE: legend.init (color) (lweight) (ltype) (label)"
+  break
+ end
+
+ local color myline mylt mylw
+ 
+ $color  = $1
+ $mylt   = $2
+ $mylw   = $3
+ $myline = $4
+
+ line -frac $legend:X $legend:Y to {$legend:X + $legend:dX} $legend:Y -c $color -lw $mylw -lt $mylt
+ textline -frac -justify 5 {$legend:X + $legend:dX + 0.01} $legend:Y "$myline"
+ $legend:Y -= $legend:dY
+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.25
+ 
+ delete -q  magBin numTotal numPoss numFound numGood
+ for mag $Mmin $Mmax $Mdel
+   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 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
+
+macro plot.camsky
+
+  dev -n sky; 
+  resize 1800 1800
+  region 194.650035  22.219991 2.5;
+  imbox o9401g0123o.stats/o9401g0123o.1808871.cm.2432793.smf
+
+  set RrefCam = Rref
+  set DrefCam = Dref
+  set Q1refCam = Q1ref
+  cplot -c blue -op 1.0 -pt cir -sz 2 RrefCam DrefCam;
+  cplot -c red  -op 0.5 -pt cir -sz 2 RrefCam DrefCam where (Q1refCam > 0.85)
+
+  data o9401g0123o.stats/o9401g0123o.mpc.fits
+  read -fits DATA Rref Dref
+  cplot -c black -pt box -sz 1 -op 0.3 Rref Dref
+
+  cplot -c blue -op 1.0 -pt box -sz 3.0 RrefCam DrefCam where (not(isnan(Xpsf)))
+  cplot -c red  -op 1.0 -pt box -sz 2.5 RrefCam DrefCam where (not(isnan(Xpsf)) && (Q1psf > 0.85))
+end
+
+############ macros used to test mpcorb_predict
+
+macro test.mpcfine
+  if ($0 != 2)
+    echo "USAGE: test.mpcfine (file)"
+    break
+  end
+
+  data $1
+  read Rr 1 Dr 2 Rc 4 Dc 5 Rp 7 Dp 8 Ra 10 Da 11 Rd 13 Dd 14
+  set dRc = 3600*(Rc - Rr)*dcos(Dr)
+  set dDc = 3600*(Dc - Dr)
+  set dRp = 3600*(Rp - Rc)*dcos(Dr)
+  set dDp = 3600*(Dp - Dc)
+  set dRa = 3600*(Ra - Rp)*dcos(Dr)
+  set dDa = 3600*(Da - Dp)
+  set dRd = 3600*(Rd - Ra)*dcos(Dr)
+  set dDd = 3600*(Dd - Da)
+
+  vstat dRc
+  vstat dRp
+  vstat dRa
+  vstat dRd
+
+end
+
+macro test.mpcorb
+  if ($0 != 3)
+    echo "USAGE: test.mpcorb (version) (options)"
+    break
+  end
+
+  local Version
+  $Version = $1
+  $Options = $2
+
+  echo mpcorb_predict moment trange.59396.v2.fits 59375.31398 {200.738197150704-2} {200.738197150704+2} {6.13117829483755-2} {6.13117829483755+2} o9375g0065o.am$Version.fits
+  exec mpcorb_predict moment trange.59396.v2.fits 59375.31398 {200.738197150704-2} {200.738197150704+2} {6.13117829483755-2} {6.13117829483755+2} o9375g0065o.am$Version.fits
+  get.camera.matches o9375g0065o.1795013.cm.2420887.smf o9375g0065o.am$Version.fits o9375g0065o.t1/o9375g0065o.t1 am$Version
+  load.camstats.glob o9375g0065o.t1/o9375g0065o.t1.XY*.am$Version.cam.stf
+
+  set dX$Version = Xpsf - Xref
+  set dY$Version = Ypsf - Yref
+
+  set dR$Version = 3600*(Rpsf - Rref)*dcos(Dref)
+  set dD$Version = 3600*(Dpsf - Dref)
+
+  plot -pt cir -sz 2 -op 0.3 dR$Version dD$Version
+end
+
+macro replot.mpcorb
+  if ($0 != 3)
+    echo "USAGE: test.mpcorb (version) (color)"
+    break
+  end
+
+  plot -pt cir -sz 2 -op 0.3 dR$1 dD$1 -c $2
+end
+
+macro load.camstats.glob.old
+ if ($0 != 2)
+   echo "USAGE: load.camstats.glob (glob)"
+   break
+ end
+ 
+ delete -q Rref Dref Mref Xref Yref Q1ref Xpsf Ypsf Rpsf Dpsf Mpsf Q1psf
+
+ list name -glob $1
+ for i 0 $name:n
+
+   data $name:$i
+   read -fits DATA -sizes
+   if ($table:Ny == 0) continue
+   read -fits DATA -list-fields -q; list tfields -join allfields; read -fits DATA $allfields
+
+   concat RrefSS Rref
+   concat DrefSS Dref
+   concat XrefSS Xref
+   concat YrefSS Yref
+   concat MrefSS Mref
+   concat Q1refSS Q1ref
+   concat X_PSF_found Xpsf
+   concat Y_PSF_found Ypsf
+   concat RA_PSF_found Rpsf
+   concat DEC_PSF_found Dpsf
+   concat M_PSF_found Mpsf
+   concat Q1_PSF_found Q1psf
+ end
+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: /branches/eam_branches/ipp-20211108/ippTests/ippCheck/configure.ac
===================================================================
--- /branches/eam_branches/ipp-20211108/ippTests/ippCheck/configure.ac	(revision 41941)
+++ /branches/eam_branches/ipp-20211108/ippTests/ippCheck/configure.ac	(revision 41941)
@@ -0,0 +1,14 @@
+AC_PREREQ(2.59)
+
+AC_INIT([ippCheck], [1.1.0], [ipp-dev@ifa.hawaii.edu])
+AC_CONFIG_SRCDIR([pantasks.pro])
+
+AM_INIT_AUTOMAKE([1.6 foreign dist-bzip2])
+AM_MAINTAINER_MODE
+
+AC_PROG_INSTALL
+
+AC_CONFIG_FILES([
+  Makefile
+])
+AC_OUTPUT
Index: /branches/eam_branches/ipp-20211108/ippTests/ippCheck/diffmask.pro
===================================================================
--- /branches/eam_branches/ipp-20211108/ippTests/ippCheck/diffmask.pro	(revision 41941)
+++ /branches/eam_branches/ipp-20211108/ippTests/ippCheck/diffmask.pro	(revision 41941)
@@ -0,0 +1,335 @@
+#!/usr/bin/env dvo
+
+$VERBOSE    = 1	       
+$COPYMASKS  = 0
+$OUTDIR     = .
+if (not($?OUTVERSION))
+  $OUTVERSION = XXX
+end
+
+$DIFF_MODE_WW = 1
+$DIFF_MODE_WS = 2
+
+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 
+$XTR_FIELDS_DIF = Xref Yref Nref Q1ref
+
+macro get.diff.eff
+ if ($0 != 5)
+  echo "USAGE: get.diff.eff (camera) (expname) (warp_id) (refsrcs.fits)"
+  echo "EX:    get.diff.eff gpc1 o9322g0253o 2367603 filename.fits"
+
+  # o9395g0075o : warp_id : 2399439, diff_id : 2156044
+  break
+ end
+
+ local dbname expname camID
+ $dbname   = $1
+ $expname  = $2
+ $warpID   = $3
+ $refsrcs  = $4
+
+  $OUTDIR = $expname.stats
+  mkdir $OUTDIR
+
+  # I have a table of sources which may be in the exposure
+  # I need to measure the mask fraction for each detection in the diffs
+  # also measure forced photometry in the diffs?  We do not actually generate a table of diff detections
+
+  # difftool -listrun to get the diff_id
+  # difftool -diffed -diff_id diffID to get the list of skyfile paths
+
+  # get the diffID for the appropriate run:
+  get.diff.run $dbname $expname $warpID $DIFF_MODE_WW
+  get.diff.skyfiles $dbname $diffID
+
+  get.diff.files
+  
+  # load the reference source table
+  # how to choose the name of the filter?
+  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"
+
+  # delete any existing output vectors
+  foreach field $SMF_FIELDS_DIF $REF_FIELDS_DIF $XTR_FIELDS_DIF
+    delete -q $field\_out
+  end
+
+  # find matching entries in cmf files and measure mask fraction
+  for i 0 $diffMask:n
+    if ("$diffCMF:$i" == "") 
+      echo "no CMF file found for $i, skipping" 
+      continue
+    end
+    get.diff.detect.stats $diffCMF:$i
+
+    if ("$diffMask:$i" == "") 
+      echo "no mask file found for $i, skipping maskstats" 
+      set Q1ref_onchp = zero(Xref_onchp) + NAN
+    else
+      get.diff.mask.stats $diffMask:$i
+    end
+
+   foreach field $SMF_FIELDS_DIF
+     concat $field\_found $field\_out
+     delete $field\_found
+   end
+   foreach field $REF_FIELDS_DIF $XTR_FIELDS_DIF
+     concat $field\_onchp $field\_out
+     delete $field\_onchp
+   end
+  end
+
+  $OUT_FIELDS = 
+  foreach field $SMF_FIELDS_DIF $REF_FIELDS_DIF $XTR_FIELDS_DIF
+    $OUT_FIELDS = $OUT_FIELDS $field\_out
+  end
+
+  sprintf output "%s/%s.%s.dif.stf" $OUTDIR $expname $OUTVERSION
+  write -fits DATA $output $OUT_FIELDS
+end
+
+macro get.diff.detect.stats
+  if ($0 != 2) 
+   echo "USAGE: get.diff.detect.stats (cmffile)"
+   echo "  assumes existence of vectors (Rref, Dref, Mref)"
+   echo "  matches known sources to detections in the cmf file"
+   echo "  calculates detection fraction as a function of magnitude"
+   echo "  also astrometric and photometric scatter vs magnitude for good detections"
+   break
+  end
+
+  local filename NAXIS1 NAXIS2
+  $filename = $1
+  if ("$filename" == "") 
+    echo "missing CMF file $filename, skipping"
+    return
+  end
+
+  # read in the image header portion
+  # rd -n SkyChip.hdr hdr $filename
+  # XXX apparently the diff cmf files have the hdr as the PHU?
+  rd hdr $filename
+  keyword hdr IMNAXIS1 NAXIS1
+  keyword hdr IMNAXIS2 NAXIS2
+
+  coords -q hdr -p       0       0; init.radec.range
+  coords -q hdr -p $NAXIS1       0; set.radec.range
+  coords -q hdr -p $NAXIS1 $NAXIS2; set.radec.range
+  coords -q hdr -p       0 $NAXIS2; set.radec.range
+
+  # generate a subset which could be on this skycell
+  # do we need to grow the range at all?
+  set keep = (Rref >= $minR) && (Rref <= $maxR) && (Dref >= $minD) && (Dref <= $maxD) 
+  foreach field $REF_FIELDS_DIF Nref
+    subset $field\_tmp = $field if keep
+  end
+
+  # convert the subset coordinates to pixels
+  set Xref_tmp = Rref_tmp
+  set Yref_tmp = Dref_tmp
+  coords -q hdr -c Xref_tmp Yref_tmp
+
+  # generate a subset of sources which are within the boundaries of the skycell
+  set keep = (Xref_tmp >= 0) && (Xref_tmp <= $NAXIS1) && (Yref_tmp >= 0) && (Yref_tmp <= $NAXIS2)
+  foreach field $REF_FIELDS_DIF Nref Xref Yref
+    subset $field\_onchp = $field\_tmp if keep
+    delete $field\_tmp
+  end
+
+  # now read the smf detections
+  data $filename
+  read -fits SkyChip.psf $SMF_FIELDS_DIF
+
+  match2d Xref_onchp Yref_onchp X_PSF Y_PSF $MATCH_RADIUS -closest
+  foreach field $SMF_FIELDS_DIF
+    reindex $field\_found = $field using index1 -keep-unmatched
+  end
+end
+
+macro get.diff.mask.stats
+  if ($0 != 2) 
+   echo "USAGE: get.diff.mask.stats (diffmask)"
+   echo "  assumes existence of vectors XrefSS YrefSS RrefSS DrefSS MrefSS" 
+   echo "  those are calculated by get.diff.detect.stats"
+   break
+  end
+
+  # NOTE: diffmask must be a real file path, not a neb path
+  local diffmask 
+  $diffmask = $1
+
+  # load the mask file
+  rd -x 0 mk $diffmask
+
+  # NOTE: assumes existence of vectors XrefSS YrefSS RrefSS DrefSS MrefSS
+  # those vectors are generated by get.diff.detect.stats
+
+  # I have a mask file, now I need to calculate something equiv to psfQf
+  set mkPerfect = (mk == 0) ? zero(mk) + 1 : zero(mk)
+
+  # I want to set an image based on the bit values but I fear
+  # the mask has been converted to a float when read
+  forcedphot mkPerfect Xref_onchp Yref_onchp -output Q1ref_onchp
+
+  subset tmp = Q1ref_onchp where (Q1ref_onchp > 0.85)
+  if ($VERBOSE) echo tmp[] of Q1ref_onchp[] are unmasked
+end
+
+## utilities
+
+macro init.radec.range
+  $minR = $RA
+  $minD = $DEC
+  $maxR = $RA
+  $maxD = $DEC
+end
+
+macro set.radec.range
+  $minR = min($minR, $RA)
+  $minD = min($minD, $DEC)
+  $maxR = max($maxR, $RA)
+  $maxD = max($maxD, $DEC)
+end
+
+macro get.diff.files
+ if ($0 != 1)
+  echo "USAGE: get.diff.files"
+  echo " uses diffPath:0 - diffPath:n"
+  echo " sets diffMask:0 - diffMask:n"
+  echo "  and diffCMF:0  - diffCMF:n"
+  break
+ end
+
+ # get the non-neb paths for the masks
+ for i 0 $diffPath:n
+   sprintf filename "%s.mask.fits" $diffPath:$i
+   break -auto off   
+   $diffMask:$i = `neb-locate -p $filename`
+   break -auto on
+   if ("$diffMask:$i" == "")
+     echo "skipping $filename"
+   else 
+     if ($COPYMASKS) exec neb-copy $filename $OUTDIR
+   end
+
+   if ($diffOrder == 0)
+     sprintf filename "%s.cmf" $diffPath:$i
+   else
+     sprintf filename "%s.inv.cmf" $diffPath:$i
+   end
+   $diffCMF:$i  = `neb-locate -p $filename`
+   if ($VERBOSE) echo "found mask: $diffMask:$i CMF: $diffCMF:$i"
+ end
+ $diffMask:n = $diffPath:n
+end
+
+macro get.diff.skyfiles
+ if ($0 != 3)
+  echo "USAGE: get.diff.skyfiles (camera) (diffID)"
+  echo "EX:    get.diff.skyfiles gpc1 2367603"
+  echo "  sets diffPath:0 - diffPath:n"
+  break
+ end
+
+ local dbname diffID
+ $dbname  = $1
+ $diffID  = $2
+
+ $command = difftool -dbname $dbname -diffskyfile -diff_id $diffID
+ 
+ queueload  myqueue -x "$command"
+ queue2book myqueue mybook
+ 
+ book npages mybook -var npages
+ 
+ if ($npages < 1) 
+  echo "ERROR: missing exposure : $command"
+  book delete mybook
+  return
+ end
+ 
+ local i
+ for i 0 $npages 
+   sprint pagename page.%03d $i
+   book getword mybook $pagename path_base -var diffPath:$i
+   book getword mybook $pagename skycell_id -var skycell:$i
+   if ($VERBOSE) echo "found diffPath $diffPath:$i"
+ end
+ $diffPath:n = $npages
+ $skycell:n = $npages
+
+ book delete mybook
+end
+
+macro get.diff.run
+ if ($0 != 5)
+  echo "USAGE: get.diff.run (camera) (expname) (warpID) (diffMode)"
+  echo "EX:    get.diff.run gpc1 o9322g0253o 2367603 1"
+  echo "  diffMode: 1 = warp-warp, 2 = warp-stack"
+  echo "  sets diffID"
+  break
+ end
+
+ local dbname expname warpID diffMode
+ $dbname   = $1
+ $expname  = $2
+ $warpID   = $3
+ $diffMode = $4
+
+ $diffID    = 0
+ $diffOrder = 0; # 0 : A of A-B, 1 : B of A-B
+
+ # we do not know if the exposure was used as A or B in the A-B diff.  First try A:
+#$command = difftool -dbname $dbname -listrun -exp_name $expname -warp_id $warpID -dist_group $distGroup -diff_mode $diffMode
+ $command = difftool -dbname $dbname -listrun -exp_name $expname -warp_id $warpID -diff_mode $diffMode
+ 
+ queueload  myqueue -x "$command"
+ queue2book myqueue mybook
+ 
+ book npages mybook -var npages
+ 
+ if ($npages < 1) 
+  # try the template versoin
+  book delete mybook
+ 
+  $command = difftool -dbname $dbname -listrun -exp_name $expname -warp_id $warpID -diff_mode $diffMode -template
+  $diffOrder = 1
+  
+  queueload  myqueue -x "$command"
+  queue2book myqueue mybook
+  
+  book npages mybook -var npages
+  
+  if ($npages < 1) 
+   echo "ERROR: missing exposure : $command"
+   book delete mybook
+   return
+  end
+ end
+ 
+ if ($npages > 1) 
+  echo "ERROR: too many exposures : $command"
+  book delete mybook
+  return
+ end
+ 
+ sprint pagename page.%03d 0
+ book getword mybook $pagename diff_id -var diffID
+ book delete mybook
+
+ if ($VERBOSE) echo "found diffID $diffID"
+end
Index: /branches/eam_branches/ipp-20211108/ippTests/ippCheck/ippCheck
===================================================================
--- /branches/eam_branches/ipp-20211108/ippTests/ippCheck/ippCheck	(revision 41941)
+++ /branches/eam_branches/ipp-20211108/ippTests/ippCheck/ippCheck	(revision 41941)
@@ -0,0 +1,146 @@
+#!/usr/bin/env dvo
+
+$PV3_CATDIR = /data/ipp105.0/eugene/3pi.dvo.masters/3pi.pv3.20170919
+
+input truth.pro
+# provides the following:
+# make.truth.sources (camera) (expname) (source) (type)
+# make.truth.sources.mpc (camera) (expname) (source)
+# make.truth.sources.dvo (camera) (expname) (source)
+# get.exposure.metadata (camera) (expname)
+
+module cammask.pro
+module warpmask.pro
+module diffmask.pro
+module camstats.pro
+module warpstats.pro
+
+## XXX I need to enable only parts of this analysis
+## for now, just comment out the warp & diff stages
+
+macro example.full.dvo
+  if ($0 != 3)
+    echo "USAGE: example.full.dvo (expname) (label)"
+    break
+  end
+
+  local expname camlabel
+  $expname = $1
+  $camlabel = $2
+ 
+  $OUTVERSION = DVO
+
+  make.truth.sources gpc1 $expname $PV3_CATDIR DVO
+  get.camera.eff     gpc1 $expname $camlabel       $REFSOURCE
+# get.warp.eff       gpc1 $expname $camID          $REFSOURCE
+# get.diff.eff       gpc1 $expname $warpID         $REFSOURCE
+
+  $KAPA = kapa -noX
+  antialias 0.5
+
+  load.camstats  $expname.stats/$expname.DVO.cam.stf
+  plot.camstats.fracs; png -name $expname.stats/$expname.dvo.fracs.png
+  plot.camstats.resid; png -name $expname.stats/$expname.dvo.resid.png
+  
+# load.warpstats $expname.stats/$expname.DVO.wrp.stf
+# png -name  $expname.stats/$expname.dvo.warphist.png
+#
+# load.camstats  $expname.stats/$expname.DVO.dif.stf
+# png -name  $expname.stats/$expname.dvo.diffhist.png
+end
+
+macro example.full.mpc
+  if ($0 != 4)
+    echo "USAGE: example.full.mpc (expname) (label) (mpcfile)"
+    break
+  end
+
+  local expname camlabel mpcfile
+  $expname = $1
+  $camlabel = $2
+  $mpcfile = $3
+ 
+  # download the most current MPCORB.DAT file with:
+  # wget ...
+
+  # generate orbit list for night / chunk with
+  # exec mpcorb_predict trange (MPCORB.DAT) (Tmin) (Tmax) (Rmin) (Rmax) (Dmin) (Dmax) (range.fits)
+  # exec mpcorb_predict trange mpc/mpcorb_extended.59396.dat 59396.3 59396.9 220.0 235.0 40.0 50.0 mpc/mpcorb.trange.59396.fits
+  # exec mpcorb_predict trange mpc/mpcorb_extended.59396.dat 59401.3 59401.9 190.0 215.0 12.0 28.0 mpc/mpcorb.trange.59401.fits
+  # visit 1 on 59401 has exposures o9401g0123o - o9401g0143o
+  # this sample has ~1900 known asteroids covering 375 sq deg. we should expect ~ 35 asteroids per exposure
+
+  # mpcorb_predict trange mpc/mpcorb_extended.59396.dat 59375.2 59375.8 190 250 -10 20 trange.59375.fits
+  # this area covers exposures o9375g0060o - o9375g0341o (visit 1) and later visits
+
+  # use a larger match radius for the asteroids
+  # 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; # sets the value of REFSOURCE 
+  get.camera.eff     gpc1 $expname $camlabel $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
+
+  $KAPA = kapa -noX
+  camstats.glob $expname.stats/$expname.*.MPC.cam.stf
+  png -name  $expname.stats/$expname.mpc.camhist.png
+  
+  warpstats.glob $expname.stats/$expname.*.MPC.warp.stf
+  png -name  $expname.stats/$expname.mpc.warphist.png
+
+  camstats.glob $expname.stats/$expname.*.MPC.diff.stf
+  png -name  $expname.stats/$expname.mpc.diffhist.png
+end
+
+macro example.old
+  make.truth.sources gpc1 o9395g0315o $PV3_CATDIR DVO
+  get.camera.eff     gpc1 o9395g0315o EU.nightlyscience test.fits
+  get.warp.eff       gpc1 o9395g0315o $camID            test.fits
+  get.diff.eff       gpc1 o9395g0315o $warpID           test.fits
+end
+
+# o9396g0312o
+# o9396g0313o
+# o9396g0314o
+# o9396g0315o
+# o9395g0315o
+
+if ($SCRIPT)
+  if ($argv:n < 3)
+    echo "USAGE: teleff.pro (mode) (expname) (label) [mpcfile]"
+    break
+  end
+
+  if ("$argv:0" == "mpc")
+    if ($argv:n != 4)
+      echo "USAGE: teleff.pro mpc (expname) (label) (mpcfile)"
+      break
+    end
+
+    example.full.mpc $argv:1 $argv:2 $argv:3
+    exit 0
+  end
+
+  if ("$argv:0" == "dvo")
+    if ($argv:n != 3)
+      echo "USAGE: teleff.pro dvo (expname) (label)"
+      break
+    end
+
+    example.full.dvo $argv:1 $argv:2
+    exit 0
+  end
+
+  echo "USAGE: teleff.pro (mode) (expname) (label) [mpcfile]"
+  echo "  mode may be 'dvo' or 'mpc'"
+  exit 2
+end
Index: /branches/eam_branches/ipp-20211108/ippTests/ippCheck/truth.pro
===================================================================
--- /branches/eam_branches/ipp-20211108/ippTests/ippCheck/truth.pro	(revision 41941)
+++ /branches/eam_branches/ipp-20211108/ippTests/ippCheck/truth.pro	(revision 41941)
@@ -0,0 +1,177 @@
+#!/usr/bin/env dvo
+
+$VERBOSE = 1	       
+$COPYMASKS = 1
+$OUTDIR = .
+
+macro make.truth.sources
+ if ($0 != 5)
+  echo "USAGE: make.truth.sources (camera) (expname) (source) (type)"
+  echo "EX:    make.truth.sources gpc1 o9322g0253o /data/ipp060.0/ipp/ippRefs/catdir.refcat.20170919.v0 DVO"
+  echo "EX:    make.truth.sources gpc1 o9322g0253o track.table.fits MPC"
+  break
+ end
+
+ if ("$4" == "MPC") 
+  make.truth.sources.mpc $1 $2 $3
+  return
+ end
+
+ if ("$4" == "DVO") 
+  make.truth.sources.dvo $1 $2 $3
+  return
+ end
+
+ echo "ERROR: unknown data source $3 ($4)"
+ break
+end
+
+# some basic parameters
+$RADIUS = 1.75; # camera radius in degrees (over-estimate)
+
+macro make.truth.sources.dvo
+ if ($0 != 4)
+  echo "USAGE: make.truth.sources.dvo (camera) (expname) (catdir)"
+  echo "EX:    make.truth.sources.dvo gpc1 o9322g0253o /data/ipp060.0/ipp/ippRefs/catdir.refcat.20170919.v0"
+  break
+ end
+
+ local dbname expname catdir
+ $dbname = $1
+ $expname = $2
+ $catdir = $3
+
+ $OUTDIR = $expname.stats
+ mkdir $OUTDIR
+
+ get.exposure.metadata $dbname $expname 
+ substr $FILTER 0 1 avecode
+ if ("$avecode" == "w")
+   $avecode = r
+ end
+
+ fprintf "%11s %10.6f %10.6f %12.6f : %s = %s" $expname $Ro $Do $MJD $FILTER $avecode
+
+ # calculate the exposure limits
+ $Dmin = $Do - $RADIUS
+ $Dmax = $Do + $RADIUS
+ 
+ if ($Do > 85)
+   echo "ERROR: exposure footprint is poorly estimated for Dec > 85 ($Do)"
+   break
+ end
+
+ $Rmin = $Ro - $RADIUS / dcos($Do)
+ $Rmax = $Ro + $RADIUS / dcos($Do)
+ 
+ catdir $catdir
+ skyregion $Rmin $Rmax $Dmin $Dmax
+
+ # -parallel needs to be set based on the catdir
+ avextract -parallel ra dec $avecode where ($avecode\:ncode > 5)
+ # avextract ra dec $avecode where (i:ncode > 5)
+ set Rref = ra
+ set Dref = dec
+ set Mref = $avecode
+
+ # the DVO truth file is specific to the exposure
+ sprintf REFSOURCE "%s/%s.DVO.fits" $OUTDIR $expname 
+ write -fits DATA $REFSOURCE Rref Dref Mref
+end
+
+macro make.truth.sources.mpc
+ if ($0 != 4)
+  echo "USAGE: make.truth.sources.mpc (camera) (expname) (mpcrange.fits)"
+  echo "EX:    make.truth.sources.mpc gpc1 o9395g0315o mpcorb.59395.fits"
+  break
+ end
+
+ local dbname expname mpcfile
+ $dbname  = $1
+ $expname = $2
+ $mpcfile = $3
+
+ $OUTDIR = $expname.stats
+ mkdir $OUTDIR
+
+ get.exposure.metadata $dbname $expname 
+ substr $FILTER 0 1 avecode
+ if ("$avecode" == "w")
+   $avecode = r
+ end
+
+ fprintf "%11s %10.6f %10.6f %12.6f : %s = %s" $expname $Ro $Do $MJD $FILTER $avecode
+
+ # for the moment, mpcorb_predict ignores precession for the first past test, 
+ # so we need to bump up the radius substantially
+ $RADIUS_BUMP = $RADIUS + 0.5
+
+ # calculate the exposure limits
+ $Dmin = $Do - $RADIUS_BUMP 
+ $Dmax = $Do + $RADIUS_BUMP
+ 
+ if ($Do > 85)
+   echo "ERROR: exposure footprint is poorly estimated for Dec > 85 ($Do)"
+   break
+ end
+
+ $Rmin = $Ro - $RADIUS_BUMP / dcos($Do)
+ $Rmax = $Ro + $RADIUS_BUMP / dcos($Do)
+ 
+ # mpcorb_predict generates a FITS table specific to the exposure with fields (Rref, Dref, Mref)
+ sprintf REFSOURCE "%s/%s.MPC.fits" $OUTDIR $expname 
+ echo mpcorb_predict moment $mpcfile $MJD $Rmin $Rmax $Dmin $Dmax $REFSOURCE
+ exec mpcorb_predict moment $mpcfile $MJD $Rmin $Rmax $Dmin $Dmax $REFSOURCE
+end
+
+macro get.exposure.metadata
+ if ($0 != 3)
+  echo "USAGE: get.exposure.metadata (camera) (expname)"
+  echo "EX:    get.exposure.metadata gpc1 o9322g0253o"
+  echo " sets Ro, Do, MJD, FILTER"
+  break
+ end
+
+ local dbname expname catdir regcommand pagename Rrad Drad
+ $dbname = $1
+ $expname = $2
+
+ $regcommand = regtool -dbname $dbname -processedexp -exp_name $expname 
+
+ queueload  regqueue -x "$regcommand"
+ queue2book regqueue regbook
+
+ $myPath = NONE
+ book npages regbook -var npages
+
+ if ($npages > 1) 
+   echo "ERROR: too many exposures for $word:0 $1 : $regcommand"
+   book delete regbook
+   next
+ end
+
+ if ($npages < 1) 
+   echo "ERROR: missing exposure for $word:0 $1 : $regcommand"
+   book delete regbook
+   next
+ end
+
+ sprint pagename page.%03d 0
+ book getword regbook $pagename ra       -var Rrad
+ book getword regbook $pagename decl     -var Drad
+ book getword regbook $pagename dateobs  -var dateraw
+ book getword regbook $pagename exp_time -var exptime
+ book getword regbook $pagename filter   -var FILTER
+
+ $TIMEREF = 2000/01/01,00:00:00
+ $TIMEFORMAT = mjd
+
+ # shift to the midpoint of the exposure
+ ctimes -abs $dateraw -var MJD
+ $MJD += 0.5*$exptime/86400.0
+
+ $Ro = $Rrad * 180 / $M_PI
+ $Do = $Drad * 180 / $M_PI
+
+ book delete regbook
+end
Index: /branches/eam_branches/ipp-20211108/ippTests/ippCheck/warpmask.pro
===================================================================
--- /branches/eam_branches/ipp-20211108/ippTests/ippCheck/warpmask.pro	(revision 41941)
+++ /branches/eam_branches/ipp-20211108/ippTests/ippCheck/warpmask.pro	(revision 41941)
@@ -0,0 +1,260 @@
+#!/usr/bin/env dvo
+
+$VERBOSE    = 1	       
+$COPYMASKS  = 0
+$OUTDIR     = .
+$OUTVERSION = XXX
+
+if (not($?USE_CAMSOURCE)) set USE_CAMSOURCE = 0
+
+$REF_FIELDS_WRP = Rref Dref Mref 
+$XTR_FIELDS_WRP = Xref Yref Nref Q1ref
+
+macro get.warp.eff
+ if ($0 != 5)
+  echo "USAGE: get.warp.eff (camera) (expname) (cam_id) (refsrcs.fits)"
+  echo "EX:    get.warp.eff gpc1 o9322g0253o 2398363 filename.fits"
+  break
+ end
+
+ local dbname expname camID
+ $dbname   = $1
+ $expname  = $2
+ $camID    = $3
+ $refsrcs  = $4
+
+  $OUTDIR = $expname.stats
+  mkdir $OUTDIR
+
+  # I have a table of sources which may be in the exposure
+  # I need to measure the mask fraction for each detection in the warps
+  # We do not actually generate a table of warp detections
+
+  # warptool -listrun to get the warp_id
+  # warptool -warped -warp_id warpID to get the list of skyfile paths
+
+  # get the warpID for the appropriate run:
+  get.warp.run $dbname $expname $camID
+  get.warp.skyfiles $dbname $warpID
+
+  get.warp.mask.files
+  
+  # load the reference source table
+  # how to choose the name of the filter?
+  # NOTE: Mref is chosen to match filter for this exposure
+  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"
+
+  # delete any existing output vectors
+  foreach field $REF_FIELDS_WRP $XTR_FIELDS_WRP
+    delete -q $field\_out
+  end
+
+  # find matching entries in smf file 
+  for i 0 $warpMask:n
+    get.warp.mask.stats $warpMask:$i
+
+    foreach field $REF_FIELDS_WRP $XTR_FIELDS_WRP
+      concat $field\_onchp $field\_out
+      delete $field\_onchp
+    end
+  end
+
+  $OUT_FIELDS = 
+  foreach field $REF_FIELDS_WRP $XTR_FIELDS_WRP
+    $OUT_FIELDS = $OUT_FIELDS $field\_out
+  end
+
+  sprintf output "%s/%s.%s.wrp.stf" $OUTDIR $expname $OUTVERSION
+  write -fits DATA $output $OUT_FIELDS
+end
+
+macro get.warp.mask.stats
+  if ($0 != 2) 
+   echo "USAGE: get.warp.mask.stats (warpmask)"
+   echo "  assumes existence of $smfreal vectors (Rref, Dref, Mref) and buffer (PHU)"
+   break
+  end
+
+  # NOTE: warpmask must be a real file path, not a neb path
+  local warpmask 
+  $warpmask = $1
+
+  # load the mask file
+  rd -x 0 mk $warpmask
+
+  # select the sources which are within the bounds of this image
+  # output is Xref_onchp Yref_onchp Rref_onchp Dref_onchp Mref_onchp
+  get.warp.detect.coords mk
+
+  # I have a mask file, now I need to calculate something equiv to psfQf
+  set mkPerfect = (mk == 0) ? zero(mk) + 1 : zero(mk)
+
+  # I want to set an image based on the bit values but I fear
+  # the mask has been converted to a float when read
+  forcedphot mkPerfect Xref_onchp Yref_onchp -output Q1ref_onchp
+
+  subset tmp = Q1ref_onchp where (Q1ref_onchp > 0.85)
+  if ($VERBOSE) echo tmp[] of Q1ref_onchp[] are unmasked
+end
+
+macro get.warp.detect.coords
+  if ($0 != 2) 
+   echo "USAGE: get.warp.detect.coords (buffer)"
+   echo "  assumes existence of $smfreal vectors (Rref, Dref, Mref) and buffer (PHU)"
+   echo "  matches known sources to detections in the smf file extensions for each chip"
+   echo "  calculates detection fraction as a function of magnitude"
+   echo "  also astrometric and photometric scatter vs magnitude for good detections"
+   break
+  end
+
+  local buffname NAXIS1 NAXIS2
+  $buffname = $1
+
+  keyword $buffname NAXIS1 NAXIS1
+  keyword $buffname NAXIS2 NAXIS2
+
+  coords -q $buffname -p       0       0; init.radec.range
+  coords -q $buffname -p $NAXIS1       0; set.radec.range
+  coords -q $buffname -p $NAXIS1 $NAXIS2; set.radec.range
+  coords -q $buffname -p       0 $NAXIS2; set.radec.range
+
+  # generate a subset which could be on this chip
+  # do we need to grow the range at all?
+  set keep = (Rref >= $minR) && (Rref <= $maxR) && (Dref >= $minD) && (Dref <= $maxD) 
+  foreach field $REF_FIELDS_WRP Nref
+    subset $field\_tmp = $field if keep
+  end
+
+  # convert the subset coordinates to pixels
+  set Xref_tmp = Rref_tmp
+  set Yref_tmp = Dref_tmp
+  coords -q $buffname -c Xref_tmp Yref_tmp
+
+  # generate a subset of sources which are within the boundaries of the chip
+  set keep = (Xref_tmp >= 0) && (Xref_tmp <= $NAXIS1) && (Yref_tmp >= 0) && (Yref_tmp <= $NAXIS2)
+  foreach field $REF_FIELDS_WRP Nref Xref Yref
+    subset $field\_onchp = $field\_tmp if keep
+    delete $field\_tmp
+  end
+end
+
+macro init.radec.range
+  $minR = $RA
+  $minD = $DEC
+  $maxR = $RA
+  $maxD = $DEC
+end
+
+macro set.radec.range
+  $minR = min($minR, $RA)
+  $minD = min($minD, $DEC)
+  $maxR = max($maxR, $RA)
+  $maxD = max($maxD, $DEC)
+end
+
+macro get.warp.mask.files
+ if ($0 != 1)
+  echo "USAGE: get.warp.mask.files"
+  echo " uses warpPath:0 - warpPath:n"
+  break
+ end
+
+ # get the non-neb paths for the masks
+ for i 0 $warpPath:n
+   sprintf filename "%s.mask.fits" $warpPath:$i
+   $warpMask:$i = `neb-locate -p $filename`
+   if ($VERBOSE) echo "found warpPath $warpMask:$i"
+   if ($COPYMASKS) exec neb-copy $filename $OUTDIR
+ end
+ $warpMask:n = $warpPath:n
+end
+
+macro get.warp.skyfiles
+ if ($0 != 3)
+  echo "USAGE: get.warp.skyfiles (camera) (warpID)"
+  echo "EX:    get.warp.skyfiles gpc1 2367603"
+  echo "  sets warpPath:0 - warpPath:n"
+  echo "  also skycell:0  - skycell:n"
+  break
+ end
+
+ local dbname warpID
+ $dbname  = $1
+ $warpID  = $2
+
+ $command = warptool -dbname $dbname -warped -warp_id $warpID
+ 
+ queueload  myqueue -x "$command"
+ queue2book myqueue mybook
+ 
+ book npages mybook -var npages
+ 
+ if ($npages < 1) 
+  echo "ERROR: missing exposure : $command"
+  book delete mybook
+  return
+ end
+ 
+ local i
+ for i 0 $npages 
+   sprint pagename page.%03d $i
+   book getword mybook $pagename path_base -var warpPath:$i
+   book getword mybook $pagename skycell_id -var skycell:$i
+   if ($VERBOSE) echo "found warpPath $warpPath:$i"
+ end
+ $warpPath:n = $npages
+ $skycell:n = $npages
+
+ book delete mybook
+end
+
+macro get.warp.run
+ if ($0 != 4)
+  echo "USAGE: get.warp.run (camera) (expname) (camID)"
+  echo "EX:    get.warp.run gpc1 o9322g0253o 2398363"
+  echo "  sets warpID"
+  break
+ end
+
+ local dbname expname camID
+ $dbname   = $1
+ $expname  = $2
+ $camID    = $3
+
+ $warpID = 0
+ $command = warptool -dbname $dbname -listrun -exp_name $expname -cam_id $camID
+ 
+ queueload  myqueue -x "$command"
+ queue2book myqueue mybook
+ 
+ book npages mybook -var npages
+ 
+ if ($npages < 1) 
+  echo "ERROR: missing exposure : $command"
+  book delete mybook
+  next
+ end
+ 
+ if ($npages > 1) 
+  echo "ERROR: too many exposures : $command"
+  book delete mybook
+  next
+ end
+ 
+ sprint pagename page.%03d 0
+ book getword mybook $pagename warp_id -var warpID
+ book delete mybook
+
+ if ($VERBOSE) echo "found warpID $warpID"
+end
Index: /branches/eam_branches/ipp-20211108/ippTests/ippCheck/warpstats.pro
===================================================================
--- /branches/eam_branches/ipp-20211108/ippTests/ippCheck/warpstats.pro	(revision 41941)
+++ /branches/eam_branches/ipp-20211108/ippTests/ippCheck/warpstats.pro	(revision 41941)
@@ -0,0 +1,107 @@
+
+macro load.warpstats
+ if ($0 != 2)
+   echo "USAGE: load.warpstats (glob)"
+   break
+ end
+ 
+ # 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
+   if ($table:Ny == 0) continue
+   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.warpstats
+   if ($i == 0)
+     set numTotal_sum = numTotal
+     set numPoss_sum  = numPoss
+   else
+     set numTotal_sum = numTotal + numTotal_sum
+     set numPoss_sum  = numPoss  + numPoss_sum
+   end
+   foreach field $FIELDS
+     concat $field\_out $field
+   end
+
+   # cleanup by deleting the loaded vectors
+   delete -q $allfields
+ end
+   
+ set numTotal = numTotal_sum
+ set numPoss  = numPoss_sum
+ delete numTotal_sum numPoss_sum
+
+ set fracPoss  = (numTotal > 0) ? numPoss  / numTotal : zero(numTotal)
+end
+
+# assumes we have loaded values from the stf tables
+# magBin, numTotal, numFound, numPoss, numGood
+macro plot.warpstats.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
+
+ # plot the detection fractions, using left axis
+ section a 0 0 1 1
+ 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 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 magBin numTotal; box $BSTY -ticks 1011 -labels 0001
+
+ plot -x hist magBin numTotal -c grey80 -lw 2
+ label +y "detected number"
+
+ vstat -q numTotal; $Ntotal = $TOTAL
+ 
+ fprintf "%d total sources" $Ntotal
+
+ legend.init 0.02 0.95 0.03 0.02
+ sprintf line "unmasked source positions (psfqf > 0.85) [%.1f pct]" {100*$FracBrightPoss};  legend.line blue  3 solid "$line"
+end
+
+# relies on the following vectors: Mref_out, Q1ref_out, PSF_QF_out, X_PSF_out
+macro calc.warpstats
+ $Mmin = 10
+ $Mmax = 25
+ $Mdel = 0.25
+ 
+ delete -q  magBin numTotal numPoss
+ for mag $Mmin $Mmax $Mdel
+   concat $mag magBin
+
+   # 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 keepMsk = keepSrc && (Q1ref_out > 0.85)
+   subset tmp = Mref_out if keepMsk
+   concat tmp[] numPoss
+ end
+
+ delete -q tmp keepSrc keepMsk keepDet
+end
