Index: /trunk/tools/eam/teleff/cammask.pro
===================================================================
--- /trunk/tools/eam/teleff/cammask.pro	(revision 41712)
+++ /trunk/tools/eam/teleff/cammask.pro	(revision 41712)
@@ -0,0 +1,342 @@
+#!/usr/bin/env dvo
+
+$VERBOSE    = 1	       
+$COPYMASKS  = 1
+$OUTDIR     = .
+$OUTVERSION = XXX
+
+if (not($?MATCH_RADIUS)) set MATCH_RADIUS = 2.0
+
+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 Rref Dref Mref
+ 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 Q1refSS = zero(RrefSS) + 1.0
+
+   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
+end
+
+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 Rref Dref Mref
+ set Nref = ramp(Rref)
+
+ # read in the PHU header
+ rd PHU $smfreal
+
+ # create the chipName list from the smffile and calculate the index numbers
+ # if a chip is missing from the smf, this may be inconsistent with chipMask:n
+ 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 
+  $chipName:$i = $word:1
+ 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
+
+   sprintf output "%s/%s.%s.%s.cam.stf" $OUTDIR $expname $chipName:$i $OUTVERSION
+   write -fits DATA $output RrefSS DrefSS MrefSS XrefSS YrefSS Q1refSS X_PSF_found Y_PSF_found M_PSF_found Q1_PSF_found Q2_PSF_found
+ end
+end
+
+macro get.chip.mask.stats
+  if ($0 != 2) 
+   echo "USAGE: get.chip.mask.stats (cammask)"
+   echo "  assumes existence of $smfreal vectors (Rref, Dref, Mref) and buffer (PHU)"
+   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 XrefSS YrefSS -output Q1refSS
+
+  subset tmp = Q1refSS where (Q1refSS > 0.85)
+
+  echo tmp[] of Q1refSS[] are unmasked
+end
+
+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
+  $chipname = $1
+  $chipnseq = $2
+
+  # read in the CHIP header
+
+  # sprintf extname "%s.hdr" $chipname
+  # rd -n $extname hdr $smfreal
+
+  rd -x $chipnseq hdr $smfreal
+
+  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) 
+  subset RrefS = Rref if keep
+  subset DrefS = Dref if keep
+  subset MrefS = Mref if keep
+
+  # convert the subset coordinates to pixels
+  set XrefS = RrefS
+  set YrefS = DrefS
+  coords -q hdr -mosaic PHU -c XrefS YrefS
+
+  # generate a subset of sources which are within the boundaries of the chip
+  set keep = (XrefS >= 0) && (XrefS <= $NAXIS1) && (YrefS >= 0) && (YrefS <= $NAXIS2)
+  subset XrefSS = XrefS if keep
+  subset YrefSS = YrefS if keep
+  subset RrefSS = RrefS if keep
+  subset DrefSS = DrefS if keep
+  subset MrefSS = MrefS if keep
+
+  # now read the smf detections
+  data $smfreal
+  #sprintf extname "%s.psf" $chipname
+  #read -fits $extname X_PSF Y_PSF RA_PSF DEC_PSF CAL_PSF_MAG PSF_QF PSF_QF_PERFECT
+
+  read -fits {$chipnseq+1} -extnum X_PSF Y_PSF RA_PSF DEC_PSF CAL_PSF_MAG PSF_QF PSF_QF_PERFECT
+
+  ## XXX note this is not very efficient because libfits needs to loop over all extensions
+  ## each time until it finds the one it wants (would be better to generate an index and 
+  ## use -fits NN -extnum
+
+  # match radius currently needs to be larger for MPC objects
+  match2d XrefSS YrefSS X_PSF Y_PSF $MATCH_RADIUS -closest
+
+  reindex X_PSF_found   = X_PSF          using index1 -keep-unmatched
+  reindex Y_PSF_found   = Y_PSF          using index1 -keep-unmatched
+  reindex M_PSF_found   = CAL_PSF_MAG    using index1 -keep-unmatched
+  reindex Q1_PSF_found  = PSF_QF         using index1 -keep-unmatched
+  reindex Q2_PSF_found  = PSF_QF_PERFECT using index1 -keep-unmatched
+  reindex RA_PSF_found  = RA_PSF         using index1 -keep-unmatched
+  reindex DEC_PSF_found = DEC_PSF        using index1 -keep-unmatched
+
+  # set dX = XrefSS - X_PSF_found; subset dXs = dX where (MrefSS < 18); vstat -q dXs
+  # set dY = YrefSS - Y_PSF_found; subset dYs = dY where (MrefSS < 18); vstat -q dYs
+  # set dM = MrefSS - M_PSF_found; subset dMs = dM where (MrefSS < 18); vstat -q dMs
+
+  # measure some stats, eg:
+  # * fraction detected vs mag, filter by PSF_QF, PSF_QF_PERFECT
+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
+
+# 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
+
+#$command = camtool -dbname $dbname -processedexp -exp_name $expname -label $camlabel
+#XXX need to generically handle the question of which version / label
+ $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
Index: /trunk/tools/eam/teleff/camstats.sh
===================================================================
--- /trunk/tools/eam/teleff/camstats.sh	(revision 41712)
+++ /trunk/tools/eam/teleff/camstats.sh	(revision 41712)
@@ -0,0 +1,259 @@
+
+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.fits 59375.31398 {200.738197150704-2} {200.738197150704+2} {6.13117829483755-2} {6.13117829483755+2} o9375g0065o.am$Version.fits $Options
+  exec mpcorb_predict moment trange.59396.fits 59375.31398 {200.738197150704-2} {200.738197150704+2} {6.13117829483755-2} {6.13117829483755+2} o9375g0065o.am$Version.fits $Options
+  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
+ if ($0 != 2)
+   echo "USAGE: load.camstats.glob (glob)"
+   break
+ end
+ 
+ delete -q Rref Dref Xref Yref Q1ref Xpsf Ypsf Rpsf Dpsf 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 Q1refSS Q1ref
+   concat X_PSF_found Xpsf
+   concat Y_PSF_found Ypsf
+   concat RA_PSF_found Rpsf
+   concat DEC_PSF_found Dpsf
+   concat Q1_PSF_found Q1psf
+ end
+end
+
+macro camstats.glob
+ if ($0 != 2)
+   echo "USAGE: camstats (glob)"
+   break
+ end
+ 
+ delete -q Rref Dref Xref Yref Q1ref Xpsf Ypsf 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
+
+   # assumes data has been loaded 
+   # results in magFound, fracFound, numFound
+   calc.camstats
+   if ($i == 0)
+     set magFound_sum = magFound
+     set numFound_sum = numFound
+     set numTotal_sum = numTotal
+     set numPoss_sum  = numPoss
+     set numOK_sum    = numOK
+   else
+     set numOK_sum    = numOK    + numOK_sum
+     set numPoss_sum  = numPoss  + numPoss_sum
+     set numFound_sum = numFound + numFound_sum
+     set numTotal_sum = numTotal + numTotal_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
+ end
+   
+ set magFound = magFound_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
+
+ threshold -q -r magFound fracFound 0.25 -range 25 {magFound[]-1}
+ $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
+ clear -s
+ resize 1800 1200
+ label -fn helvetica 18
+ 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
+ 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
+
+ section b 0 0 1 1
+ lim magFound numFound; box $BSTY -ticks 1011 -labels 0001; plot -x hist magFound numFound -c grey70 -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
+
+ $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
+ line -frac $Xlegend $Ylegend to {$Xlegend + $dXlegend} $Ylegend -c red70 -lw 3; textline -frac -justify 5 {$Xlegend + $dXlegend + 0.01} $Ylegend "all sources detected"                    ; $Ylegend -= $dYlegend
+ line -frac $Xlegend $Ylegend to {$Xlegend + $dXlegend} $Ylegend -c red   -lw 3; textline -frac -justify 5 {$Xlegend + $dXlegend + 0.01} $Ylegend "sources detected with psfqf > 0.85"
+end
+
+macro calc.camstats
+ $Mmin = 10
+ $Mmax = 25
+ $Mdel = 0.25
+ 
+ delete -q fracFound magFound numFound numTotal numPoss fracPoss numOK fracOK
+ 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 $NPTS numFound
+
+   # 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
+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
Index: /trunk/tools/eam/teleff/diffmask.pro
===================================================================
--- /trunk/tools/eam/teleff/diffmask.pro	(revision 41712)
+++ /trunk/tools/eam/teleff/diffmask.pro	(revision 41712)
@@ -0,0 +1,326 @@
+#!/usr/bin/env dvo
+
+$VERBOSE    = 1	       
+$COPYMASKS  = 1
+$OUTDIR     = .
+$OUTVERSION = XXX
+
+$DIFF_MODE_WW = 1
+$DIFF_MODE_WS = 2
+
+if (not($?MATCH_RADIUS)) set MATCH_RADIUS = 2.0
+
+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; read -fits DATA Rref Dref Mref
+  set Nref = ramp(Rref)
+  if ($VERBOSE) echo "loaded Rref[] sources from reference $refsrcs"
+
+  # 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 Q1refSS = zero(XrefSS) + NAN
+    else
+      get.diff.mask.stats $diffMask:$i
+    end
+
+    sprintf output "%s/%s.%s.%s.diff.stf" $OUTDIR $expname $skycell:$i $OUTVERSION
+    write -fits DATA $output RrefSS DrefSS MrefSS XrefSS YrefSS Q1refSS X_PSF_found Y_PSF_found M_PSF_found Q1_PSF_found Q2_PSF_found
+  end
+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
+  $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) 
+  subset RrefS = Rref if keep
+  subset DrefS = Dref if keep
+  subset MrefS = Mref if keep
+
+  # convert the subset coordinates to pixels
+  set XrefS = RrefS
+  set YrefS = DrefS
+  coords -q hdr -c XrefS YrefS
+
+  # generate a subset of sources which are within the boundaries of the skycell
+  set keep = (XrefS >= 0) && (XrefS <= $NAXIS1) && (YrefS >= 0) && (YrefS <= $NAXIS2)
+  subset XrefSS = XrefS if keep
+  subset YrefSS = YrefS if keep
+  subset RrefSS = RrefS if keep
+  subset DrefSS = DrefS if keep
+  subset MrefSS = MrefS if keep
+
+  # now read the smf detections
+  data $filename
+  read -fits SkyChip.psf X_PSF Y_PSF RA_PSF DEC_PSF CAL_PSF_MAG PSF_QF PSF_QF_PERFECT
+
+  match2d XrefSS YrefSS X_PSF Y_PSF $MATCH_RADIUS -closest
+
+  reindex X_PSF_found = X_PSF using index1 -keep-unmatched
+  reindex Y_PSF_found = Y_PSF using index1 -keep-unmatched
+  reindex M_PSF_found = CAL_PSF_MAG using index1 -keep-unmatched
+  reindex Q1_PSF_found = PSF_QF using index1 -keep-unmatched
+  reindex Q2_PSF_found = PSF_QF_PERFECT using index1 -keep-unmatched
+
+  set dX = XrefSS - X_PSF_found; subset dXs = dX where (MrefSS < 18); vstat -q dXs
+  set dY = YrefSS - Y_PSF_found; subset dYs = dY where (MrefSS < 18); vstat -q dYs
+  set dM = MrefSS - M_PSF_found; subset dMs = dM where (MrefSS < 18); vstat -q dMs
+
+  # measure some stats, eg:
+  # * fraction detected vs mag, filter by PSF_QF, PSF_QF_PERFECT
+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 XrefSS YrefSS -output Q1refSS
+
+  subset tmp = Q1refSS where (Q1refSS > 0.85)
+
+  if ($VERBOSE) echo tmp[] of Q1refSS[] are unmasked
+end
+
+macro get.diff.detect.coords
+  if ($0 != 2) 
+   echo "USAGE: get.diff.detect.coords (buffer)"
+   echo "  assumes existence of 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
+  $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) 
+  subset RrefS = Rref if keep
+  subset DrefS = Dref if keep
+  subset MrefS = Mref if keep
+
+  # convert the subset coordinates to pixels
+  set XrefS = RrefS
+  set YrefS = DrefS
+  coords -q $buffname -c XrefS YrefS
+
+  # generate a subset of sources which are within the boundaries of the chip
+  set keep = (XrefS >= 0) && (XrefS <= $NAXIS1) && (YrefS >= 0) && (YrefS <= $NAXIS2)
+  subset XrefSS = XrefS if keep
+  subset YrefSS = YrefS if keep
+  subset RrefSS = RrefS if keep
+  subset DrefSS = DrefS if keep
+  subset MrefSS = MrefS if keep
+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.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
+   $diffMask:$i = `neb-locate -p $filename`
+   if ($COPYMASKS) exec neb-copy $filename $OUTDIR
+
+   sprintf filename "%s.cmf" $diffPath:$i
+   $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
+#$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) 
+  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 diff_id -var diffID
+ book delete mybook
+
+ if ($VERBOSE) echo "found diffID $diffID"
+end
Index: /trunk/tools/eam/teleff/teleff.pro
===================================================================
--- /trunk/tools/eam/teleff/teleff.pro	(revision 41712)
+++ /trunk/tools/eam/teleff/teleff.pro	(revision 41712)
@@ -0,0 +1,129 @@
+#!/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)
+
+input cammask.pro
+input warpmask.pro
+input diffmask.pro
+
+input camstats.sh
+input warpstats.sh
+
+macro example.full.dvo
+  if ($0 != 2)
+    echo "USAGE: example.full.dvo (expname)"
+    break
+  end
+
+  local expname
+  $expname = $1
+ 
+  $OUTVERSION = DVO
+
+  make.truth.sources gpc1 $expname $PV3_CATDIR DVO
+  get.camera.eff     gpc1 $expname EU.nightlyscience $REFSOURCE
+  get.warp.eff       gpc1 $expname $camID            $REFSOURCE
+  get.diff.eff       gpc1 $expname $warpID           $REFSOURCE
+
+  $KAPA = kapa -noX
+  camstats.glob $expname.stats/$expname.*.DVO.cam.stf
+  png -name  $expname.stats/$expname.dvo.camhist.png
+  
+  warpstats.glob $expname.stats/$expname.*.DVO.warp.stf
+  png -name  $expname.stats/$expname.dvo.warphist.png
+
+  camstats.glob $expname.stats/$expname.*.DVO.diff.stf
+  png -name  $expname.stats/$expname.dvo.diffhist.png
+end
+
+macro example.full.mpc
+  if ($0 != 3)
+    echo "USAGE: example.full.mpc (expname) (mpcfile)"
+    break
+  end
+
+  local expname mpcfile
+  $expname = $1
+  $mpcfile = $2
+ 
+  # 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
+
+  # use a larger match radius for the asteroids
+  $MATCH_RADIUS = 3.0
+  $OUTVERSION = MPC
+  $COPYMASKS = 0
+
+  make.truth.sources gpc1 $expname $mpcfile MPC
+  get.camera.eff     gpc1 $expname EU.nightlyscience $REFSOURCE
+  get.warp.eff       gpc1 $expname $camID            $REFSOURCE
+  get.diff.eff       gpc1 $expname $warpID           $REFSOURCE
+
+  $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 < 2)
+    echo "USAGE: teleff.pro (mode) (expname) [mpcfile]"
+    break
+  end
+
+  if ("$argv:0" == "mpc")
+    if ($argv:n != 3)
+      echo "USAGE: teleff.pro mpc (expname) (mpcfile)"
+      break
+    end
+
+    example.full.mpc $argv:1 $argv:2
+    exit 0
+  end
+
+  if ("$argv:0" == "dvo")
+    if ($argv:n != 2)
+      echo "USAGE: teleff.pro dvo (expname)"
+      break
+    end
+
+    example.full.dvo $argv:1
+    exit 0
+  end
+
+  echo "USAGE: teleff.pro (mode) (expname) [mpcfile]"
+  echo "  mode may be 'dvo' or 'mpc'"
+  exit 2
+end
Index: /trunk/tools/eam/teleff/truth.pro
===================================================================
--- /trunk/tools/eam/teleff/truth.pro	(revision 41712)
+++ /trunk/tools/eam/teleff/truth.pro	(revision 41712)
@@ -0,0 +1,174 @@
+#!/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) (subset.fits)"
+  echo "EX:    make.truth.sources.mpc gpc1 o9395g0315o subset.chunk.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 filter  -var FILTER
+
+ $TIMEREF = 2000/01/01,00:00:00
+ $TIMEFORMAT = mjd
+
+ ctimes -abs $dateraw -var MJD
+
+ $Ro = $Rrad * 180 / $M_PI
+ $Do = $Drad * 180 / $M_PI
+
+ book delete regbook
+end
Index: /trunk/tools/eam/teleff/warpmask.pro
===================================================================
--- /trunk/tools/eam/teleff/warpmask.pro	(revision 41712)
+++ /trunk/tools/eam/teleff/warpmask.pro	(revision 41712)
@@ -0,0 +1,234 @@
+#!/usr/bin/env dvo
+
+$VERBOSE    = 1	       
+$COPYMASKS  = 1
+$OUTDIR     = .
+$OUTVERSION = XXX
+
+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?
+  data $refsrcs
+  # NOTE: Mref is chosen to match filter for this exposure
+  data $refsrcs; read -fits DATA Rref Dref Mref
+  set Nref = ramp(Rref)
+  if ($VERBOSE) echo "loaded Rref[] sources from reference $refsrcs"
+
+  # find matching entries in smf file 
+  for i 0 $warpMask:n
+    get.warp.mask.stats $warpMask:$i
+
+    sprintf output "%s/%s.%s.%s.warp.stf" $OUTDIR $expname $skycell:$i $OUTVERSION
+    write -fits DATA $output RrefSS DrefSS MrefSS XrefSS YrefSS Q1refSS
+  end
+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 XrefSS YrefSS RrefSS DrefSS MrefSS
+  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 XrefSS YrefSS -output Q1refSS
+
+  subset tmp = Q1refSS where (Q1refSS > 0.85)
+
+  if ($VERBOSE) echo tmp[] of Q1refSS[] 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
+  $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) 
+  subset RrefS = Rref if keep
+  subset DrefS = Dref if keep
+  subset MrefS = Mref if keep
+
+  # convert the subset coordinates to pixels
+  set XrefS = RrefS
+  set YrefS = DrefS
+  coords -q $buffname -c XrefS YrefS
+
+  # generate a subset of sources which are within the boundaries of the chip
+  set keep = (XrefS >= 0) && (XrefS <= $NAXIS1) && (YrefS >= 0) && (YrefS <= $NAXIS2)
+  subset XrefSS = XrefS if keep
+  subset YrefSS = YrefS if keep
+  subset RrefSS = RrefS if keep
+  subset DrefSS = DrefS if keep
+  subset MrefSS = MrefS if keep
+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: /trunk/tools/eam/teleff/warpstats.sh
===================================================================
--- /trunk/tools/eam/teleff/warpstats.sh	(revision 41712)
+++ /trunk/tools/eam/teleff/warpstats.sh	(revision 41712)
@@ -0,0 +1,119 @@
+
+macro warpstats.glob
+ if ($0 != 2)
+   echo "USAGE: warpstats (glob)"
+   break
+ end
+ 
+ delete -q Rref Dref Xref Yref Q1ref
+
+ 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
+
+   # assumes data has been loaded 
+   # results in magFound, fracFound, numFound
+   calc.warpstats
+   if ($i == 0)
+     set magFound_sum = magFound
+     set numTotal_sum = numTotal
+     set numPoss_sum  = numPoss
+   else
+     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
+ end
+   
+ set magFound = magFound_sum
+ set numTotal = numTotal_sum
+ set fracPoss  = numPoss_sum  / numTotal_sum
+
+ plot.warpstats
+
+ # plot.warpsky
+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 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
+ 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
+ line -c blue -lt dot -lw 2 100 $FracBrightPoss to 0 $FracBrightPoss
+ fprintf "%.1f-pct of sources are unmasked" {100*$FracBrightPoss}
+
+ 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"
+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
+
+
+macro calc.warpstats
+ $Mmin = 10
+ $Mmax = 25
+ $Mdel = 0.25
+ 
+ delete -q magFound numTotal numPoss fracPoss
+ for mag $Mmin $Mmax $Mdel
+   concat $mag magFound
+
+   # 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
+
+   # 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
+ end
+end
+