Index: /trunk/tools/eam/rawfix.20230221/src/getvidstats.sh
===================================================================
--- /trunk/tools/eam/rawfix.20230221/src/getvidstats.sh	(revision 42743)
+++ /trunk/tools/eam/rawfix.20230221/src/getvidstats.sh	(revision 42743)
@@ -0,0 +1,304 @@
+#!/usr/bin/env mana
+$EXT = vdx
+$OUTDIR = vdx.v2
+$KEEP_OLD = 0
+
+# get.vidstats.expname (expname)
+macro get.vidstats.expname
+  if ($0 != 2)
+    echo "USAGE: get.vidstats.expname (expname)"
+    break
+  end
+  
+  local expname outfile
+  $expname = $1
+  
+  mkdir $OUTDIR
+  $outfile = $OUTDIR/$expname.$EXT
+
+  # given expname, get list of uri from rawExp/rawImfile
+  queueload  regqueue -x "regtool -dbname gpc1 -processedimfile -exp_name $expname"
+  queue2book regqueue regbook
+  
+  book npages regbook -var npages
+  if ($npages == 0) 
+    echo "ERROR: no chip imfiles for $expname"
+    break
+  end
+  
+  if ($npages != 60) 
+    echo "WARNING: unexpected number of chip imfiles: $npages vs 60"
+  end
+  
+  if ($KEEP_OLD)
+    file $outfile isFound
+    if ($isFound) return
+  else
+    exec rm -f $outfile
+  end
+
+  # define the header
+  name.values
+  $HAS_HEADER = 0
+
+  local i pagename uri
+  for i 0 $npages
+    sprint pagename page.%03d $i
+    book getword regbook $pagename uri -var uri
+  
+    ## NOTE : returned variable 'filename' must be global 
+    get.realfile.from.nebfile $uri filename
+    if ($EXECSTATUS)
+      echo "missing $uri, skipping"
+      continue
+    end
+
+    echo "trying $filename"
+    get.vidstats.imfile $filename
+
+    ## HAS_VIDEO is set by get.vidstats.imfile
+    if ($HAS_VIDEO)
+      book getword regbook $pagename class_id -var chipname
+
+      output $outfile
+      # only write out the header once, first line in the file 
+      if (not($HAS_HEADER))
+        echo "# chipname $VIDNAMES"  
+	$HAS_HEADER = 1
+      end
+      echo "$chipname $VIDSTATS"
+      output stdout
+    end
+  end
+  book delete regbook
+end
+
+macro get.vidstats.imfile
+ if ($0 != 2)
+  echo "USAGE: get.vidstats.imfile (file)"
+  break
+ end
+
+ local filename
+ $filename = $1
+
+ init.vidcells
+ $HAS_VIDEO = 0
+
+ # to find the video cell, we can look at the PHU and grab CELLMODE
+ rd hdr $filename
+ keyword hdr CELLMODE cellmode
+ strstr -q $cellmode V -start Nstart -length Nchar
+ if ($Nstart == -1)
+   # echo "no video data for this chip"
+   return
+ end
+
+ sprintf extname "video_table_xy%02d" cellXY[$Nstart]
+
+ # read the video cell table (need cnpix1,2 to measure flat-corrected flux)
+ data $filename
+ break -auto off
+ read -fits $extname frame_start_time frame_complete_time centroid_x centroid_y fwhm_x fwhm_y flux cnpix1 cnpix2
+ if (not($STATUS))
+   echo "missing FITS table field for $filename"
+   return
+ end
+ break -auto on
+
+ # set default values for outputs
+ init.values
+
+ # calculate the time offset
+ set dtime = frame_start_time - frame_start_time[0]
+ set dx = centroid_x - cnpix1
+ set dy = centroid_y - cnpix2
+
+ # a good star measurement is within the window 15 < (dx,dy) 85 pixels .
+
+ set good = (dtime > 4) && (dx > 15) && (dx < 85) && (dy > 15) && (dy < 85) && (dtime < dtime[-1] - 2.0)
+ subset dts    = dtime where good
+ subset dxs    = dx    where good
+ subset dys    = dy    where good
+ subset flux_s = flux  where good
+ $GOODFRAC = dxs[] / dx[]
+
+ if (dxs[] == 0)
+   save.values
+   return
+ end
+
+ # generate percentiles: (5, 16, 50, 84, 95) for dxs, dys
+ get.percentiles dxs; $DX_SIG_RAW = $dPCT_68 / 2.0; $DX_D90_RAW = $dPCT_90; $DX_MED_RAW = $PCT_50
+ get.percentiles dys; $DY_SIG_RAW = $dPCT_68 / 2.0; $DY_D90_RAW = $dPCT_90; $DY_MED_RAW = $PCT_50
+ vstat -q flux_s; $FLUX_MEDIAN = $MEDIAN; $FLUX_STDEV = $SIGMA
+
+ # smoothed version of dxs, dys:
+ vmedfilt dxs dxsmooth 10
+ vmedfilt dys dysmooth 10
+ set seqsm = ramp(dxsmooth)
+
+ # smoothed vector stdev
+ get.percentiles dxsmooth; $DX_SIG_SMT = $dPCT_68 / 2.0; $DX_D90_SMT = $dPCT_90; $DX_MED_SMT = $PCT_50
+ get.percentiles dysmooth; $DY_SIG_SMT = $dPCT_68 / 2.0; $DY_D90_SMT = $dPCT_90; $DY_MED_SMT = $PCT_50
+
+ # measure long-term drift
+ vstat -q dxsmooth where (seqsm < 10);             $DX_DRIFT_S = $MEDIAN
+ vstat -q dxsmooth where (seqsm > seqsm[-1] - 10); $DX_DRIFT_E = $MEDIAN
+ $DX_DRIFT = $DX_DRIFT_E - $DX_DRIFT_S
+
+ vstat -q dysmooth where (seqsm < 10);             $DY_DRIFT_S = $MEDIAN
+ vstat -q dysmooth where (seqsm > seqsm[-1] - 10); $DY_DRIFT_E = $MEDIAN
+ $DY_DRIFT = $DY_DRIFT_E - $DY_DRIFT_S
+ 
+ # scatter out of min range
+ vstat -q dxsmooth where (abs(dxsmooth - $DX_MED_SMT) > 2.0); $FX_OUT_SMT = $NPIX / dxsmooth[]
+ vstat -q dysmooth where (abs(dysmooth - $DY_MED_SMT) > 2.0); $FY_OUT_SMT = $NPIX / dysmooth[]
+
+ save.values
+end
+
+#### Utilities ####
+
+macro init.values
+ $GOODFRAC    = 0
+ $FLUX_MEDIAN = NAN
+ $FLUX_STDEV  = NAN
+ $DX_SIG_RAW  = NAN
+ $DY_SIG_RAW  = NAN
+ $DX_D90_RAW  = NAN
+ $DY_D90_RAW  = NAN
+ $DX_SIG_SMT  = NAN
+ $DY_SIG_SMT  = NAN
+ $DX_D90_SMT  = NAN
+ $DY_D90_SMT  = NAN
+ $DX_DRIFT    = NAN
+ $DY_DRIFT    = NAN
+ $FX_OUT_SMT  = NAN
+ $FY_OUT_SMT  = NAN
+end
+
+# output values are all Global
+macro save.values
+ $HAS_VIDEO = 1
+
+ sprintf VIDSTAT1 "%5.3f %9.1f %9.1f" $GOODFRAC $FLUX_MEDIAN $FLUX_STDEV
+ sprintf VIDSTAT2 "%6.2f %6.2f : %6.2f %6.2f" $DX_SIG_RAW $DY_SIG_RAW $DX_D90_RAW $DY_D90_RAW 
+ sprintf VIDSTAT3 "%6.2f %6.2f : %6.2f %6.2f" $DX_SIG_SMT $DY_SIG_SMT $DX_D90_SMT $DY_D90_SMT
+ sprintf VIDSTAT4 "%6.2f %6.2f : %6.4f %6.4f" $DX_DRIFT $DY_DRIFT $FX_OUT_SMT $FY_OUT_SMT
+ sprintf VIDSTATS "%s : %s : %s : %s" "$VIDSTAT1" "$VIDSTAT2" "$VIDSTAT3" "$VIDSTAT4"
+end
+
+# output values are all Global
+macro name.values
+ sprintf VIDNAME1 "%5s %9s %9s" GOODFRAC FLUX_MEDIAN FLUX_STDEV
+ sprintf VIDNAME2 "%6s %6s : %6s %6s" DX_SIG_RAW DY_SIG_RAW DX_D90_RAW DY_D90_RAW 
+ sprintf VIDNAME3 "%6s %6s : %6s %6s" DX_SIG_SMT DY_SIG_SMT DX_D90_SMT DY_D90_SMT
+ sprintf VIDNAME4 "%6s %6s : %5s %5s" DX_DRIFT DY_DRIFT FX_OUT_SMT FY_OUT_SMT
+ sprintf VIDNAMES "%s : %s : %s : %s" "$VIDNAME1" "$VIDNAME2" "$VIDNAME3" "$VIDNAME4"
+end
+
+# only used if needed
+macro init.vidcells
+  if (cellXY?[]) return
+
+  local ix iy cellXY
+
+  # create a lookup table sequence -> xy
+  delete -q cellXY
+  for iy 0 8
+    for ix 0 8
+      $cellXY = $ix\$iy
+      concat $cellXY cellXY
+    end
+  end
+end
+
+# set minRange - maxRange as 2 x (0.05 to 0.95) range
+macro get.threshrange
+  if ($0 != 2)
+    echo "USAGE: get.threshrange (vector)"
+    break
+  end
+
+  local vecname
+  $vecname = $1
+
+  vstat -q $vecname
+  histogram $vecname _tmpNx 0 {2*$MEDIAN} {0.01*$MEDIAN} -range _tmpdx
+  cumulative _tmpNx _tmpSx 
+  set _tmpFx = _tmpSx / _tmpSx[-1]
+  threshold -saturate -q _tmpdx _tmpFx 0.05 
+  $minT = $threshold
+
+  threshold -saturate -q _tmpdx _tmpFx 0.95
+  $maxT = $threshold
+  
+  $allRange = $maxT - $minT
+  $midRange = 0.5*($maxT + $minT)
+  $minRange = $midRange - $allRange
+  $maxRange = $midRange + $allRange
+end
+
+# set minRange - maxRange as 2 x (0.05 to 0.95) range
+macro get.percentiles
+  if ($0 != 2)
+    echo "USAGE: get.percentiles (vector)"
+    break
+  end
+
+  local vecname
+  $vecname = $1
+
+  # percentiles to calculate:
+  # 16, 84 (1-sigma)
+  #  5, 95 (2-sigma)
+  # 50 (median)
+  vstat -q $vecname
+  if ($MIN == $MAX)
+    $PCT_05 = $MIN
+    $PCT_16 = $MIN
+    $PCT_50 = $MIN
+    $PCT_84 = $MIN
+    $PCT_95 = $MIN
+  else
+
+    histogram $vecname _tmpNx $MIN $MAX {($MAX - $MIN) / 1000} -range _tmpdx
+    cumulative _tmpNx _tmpSx 
+    set _tmpFx = _tmpSx / _tmpSx[-1]
+    
+    threshold -saturate -q _tmpdx _tmpFx 0.05; $PCT_05 = $threshold
+    threshold -saturate -q _tmpdx _tmpFx 0.16; $PCT_16 = $threshold
+    threshold -saturate -q _tmpdx _tmpFx 0.50; $PCT_50 = $threshold
+    threshold -saturate -q _tmpdx _tmpFx 0.84; $PCT_84 = $threshold
+    threshold -saturate -q _tmpdx _tmpFx 0.95; $PCT_95 = $threshold
+  end
+    
+  $dPCT_90 = $PCT_95 - $PCT_05
+  $dPCT_68 = $PCT_84 - $PCT_16
+end
+
+## NOTE : output filename must be global 
+macro get.realfile.from.nebfile
+  if ($0 != 3)
+    echo "USAGE: get.realfile.from.nebfile (nebname) (realname)"
+    break
+  end
+
+  $$2 = `neb-locate -p $1`
+  if ($EXECSTATUS > 0)
+    $$2 = NONE
+  end
+end
+
+#### run as script
+if ($SCRIPT)
+  if ($argv:n != 1)
+    echo "USAGE: getvidstats.sh (expname) : output to $OUTDIR/expname.$EXT"
+    exit 2
+  end
+
+  get.vidstats.expname $argv:0
+  exit 0
+end
+