Changeset 36677

Timestamp:

Apr 20, 2014, 5:52:37 PM (12 years ago)

Author:

eugene

Message:

merge changes from trunk

Location:

branches/eam_branches/ipp-20140402

Files:

• . (modified) (1 prop)
• Ohana/src/libohana/src (modified) (1 prop)
• PS-IPP-PStamp/lib/PS/IPP/PStamp/Job.pm (modified) (2 diffs)
• ippScripts/scripts/dist_defineruns.pl (modified) (4 diffs)
• ippScripts/scripts/ipp_apply_burntool_single.pl (modified) (1 prop)
• ippScripts/scripts/psphot_fullforce_summary.pl (modified) (1 diff)
• ippScripts/scripts/publish_file.pl (modified) (1 prop)
• ippTasks/pstamp.pro (modified) (1 diff)
• ippTools/share/disttool_definebyquery_SSdiff.sql (modified) (1 diff)
• ippTools/share/disttool_definebyquery_camera.sql (modified) (1 diff)
• ippTools/share/disttool_definebyquery_chip_bg.sql (modified) (1 diff)
• ippTools/share/disttool_definebyquery_diff.sql (modified) (1 diff)
• ippTools/share/disttool_definebyquery_fake.sql (modified) (1 diff)
• ippTools/share/disttool_definebyquery_ff.sql (modified) (1 diff)
• ippTools/share/disttool_definebyquery_raw.sql (modified) (1 diff)
• ippTools/share/disttool_definebyquery_sky.sql (modified) (1 diff)
• ippTools/share/disttool_definebyquery_sky_singlefilter.sql (modified) (1 diff)
• ippTools/share/disttool_definebyquery_skycal.sql (modified) (2 diffs)
• ippTools/share/disttool_definebyquery_stack.sql (modified) (1 diff)
• ippTools/share/disttool_definebyquery_warp.sql (modified) (1 diff)
• ippTools/share/disttool_definebyquery_warp_bg.sql (modified) (1 diff)
• ippTools/share/warptool_warped.sql (modified) (1 diff)
• ippTools/src/fftool.c (modified) (1 diff)
• ippTools/src/fftoolConfig.c (modified) (1 diff)
• ippconfig (modified) (1 prop)
• ippconfig/gpc1/ppStack.config (modified) (2 diffs)
• ippconfig/recipes/nightly_science.config (modified) (2 diffs)
• ippconfig/recipes/ppStack.config (modified) (1 diff)
• ippconfig/recipes/ppSub.config (modified) (1 prop)
• ippconfig/recipes/psphot.config (modified) (1 diff)
• ippconfig/recipes/reductionClasses.mdc (modified) (1 diff, 1 prop)
• ppBackground/src/ppBackgroundStack.h (modified) (1 diff)
• ppBackground/src/ppBackgroundStackCamera.c (modified) (4 diffs)
• ppBackground/src/ppBackgroundStackData.c (modified) (1 diff)
• ppBackground/src/ppBackgroundStackLoop.c (modified) (8 diffs)
• ppBackground/src/ppBackgroundStackMath.c (modified) (8 diffs)
• ppSkycell/src/ppSkycellLoop.c (modified) (1 diff)
• psModules/src/imcombine/pmStack.c (modified) (13 diffs)
• psModules/src/imcombine/pmSubtraction.c (modified) (1 diff)
• psModules/src/objects/pmPSFtryFitPSF.c (modified) (3 diffs)
• psastro/src/psastroLoadCrosstalk.c (modified) (8 diffs)
• psphot/src (modified) (1 prop)
• psphot/src/psphotFullForceReadout.c (modified) (2 diffs)
• psphot/src/psphotFullForceSummary.c (modified) (1 diff)
• psphot/src/psphotLensing.c (modified) (1 diff)
• psphot/src/psphotStackImageLoop.c (modified) (1 prop)
• pstamp/scripts/pstamp_save_server_status.pl (modified) (1 prop)
• tools/dvodb_calib/README (modified) (1 diff)
• tools/dvodb_calib/gene.run.relphot.3pinew (copied) (copied from trunk/tools/dvodb_calib/gene.run.relphot.3pinew )
• tools/dvodb_calib/gene.run.relphot.mdnew (copied) (copied from trunk/tools/dvodb_calib/gene.run.relphot.mdnew )

branches/eam_branches/ipp-20140402/PS-IPP-PStamp/lib/PS/IPP/PStamp/Job.pm

@@ -468 +468 @@
             $skycell_id = $component;
             $choose_components = 1;
+            if ($req_type eq 'byskycell' and $skycenter) {
+                # Add a wide coordinate cut to the releasetool command.
+                # This greatly speeds up the query by reducing the number of exposures considered
+                # to those near the region of interest.
+                # The actual radius value that we use is not critical because the skycell cut does the relevant
+                # selection.
+                $command .= " -ra $row->{CENTER_X} -decl $row->{CENTER_Y} -radius 5.0";
+            }
         } else {
             $command = "$warptool -pstamp_order -dbname $imagedb";
@@ -1094 +1102 @@
             setErrorCodes($rowList, $PSTAMP_NO_IMAGE_MATCH);
         } else {
-            # XXX: We are not applying the fwhm cuts
+            # XXX: We are not applying the fwhm cuts unless we use releasetool
             foreach my $skycell (@$skycells) {
                 my $these_results = lookup($ipprc, $rowList, $imagedb, "byskycell", $stage, undef,

branches/eam_branches/ipp-20140402/ippScripts/scripts/dist_defineruns.pl

@@ -42 +42 @@
 
 # Parse the command-line arguments
-my ($stage, $stage_limit, $dist_root, $workdir, $no_magic);
+my ($stage, $stage_limit, $dist_root, $workdir, $no_magic, $rerun);
 my ($dbname, $save_temps, $verbose, $no_update, $logfile);
 my @labels;
@@ -51 +51 @@
            'workdir=s'      => \$workdir,     # workdir
            'stage_limit=s'  => \$stage_limit,# maximum number of runs queued for each stage
+           'rerun'          => \$rerun,      # queue new runs even if one exists
 #           'dist_root=s'    => \$dist_root,  # root of distribution work area
            'no_magic'       => \$no_magic,   # queue runs without requiring magic (for testing only)
@@ -56 +57 @@
            'verbose'        => \$verbose,    # Print stuff?
            'no-update'      => \$no_update,  # Don't update the database
+           'pretend'        => \$no_update,  # Don't update the database
            'save-temps'     => \$save_temps, # Save temporary files?
            'logfile=s'      => \$logfile,
@@ -104 +106 @@
         $command .= " -singlefilter" if $single;
         $command .= " -no_magic" if $no_magic;
+        $command .= " -rerun" if $rerun;
         $command .= " -pretend" if $no_update;
         $command .= " -limit $stage_limit" if $stage_limit;

branches/eam_branches/ipp-20140402/ippScripts/scripts/psphot_fullforce_summary.pl

@@ -85 +85 @@
     my $results;
     {
-        my $command = "$fftool -result -ff_id $ff_id";
+        my $command = "$fftool -result -ff_id $ff_id -quality 0";
         $command .= " -dbname $dbname" if defined $dbname;
         my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =

branches/eam_branches/ipp-20140402/ippTasks/pstamp.pro

@@ -951 +951 @@
     # overloading nfs
 
-    npending    10
+    npending    5
 
     task.exec

branches/eam_branches/ipp-20140402/ippTools/share/disttool_definebyquery_SSdiff.sql

@@ -22 +22 @@
       AND distTarget.state = 'enabled'
       AND rcInterest.state = 'enabled'
-      AND distRun.dist_id IS NULL
-
       AND ((diffRun.state = 'full') OR (distTarget.clean AND diffRun.state = 'cleaned'))

branches/eam_branches/ipp-20140402/ippTools/share/disttool_definebyquery_camera.sql

@@ -21 +21 @@
 WHERE distTarget.state = 'enabled'
     AND rcInterest.state = 'enabled'
-    AND distRun.dist_id IS NULL      -- no existing distRun 
     AND ((camRun.state = 'full') OR (distTarget.clean AND camRun.state = 'cleaned'))

branches/eam_branches/ipp-20140402/ippTools/share/disttool_definebyquery_chip_bg.sql

@@ -23 +23 @@
 WHERE distTarget.state = 'enabled'
     AND rcInterest.state = 'enabled'
-    AND distRun.dist_id IS NULL      -- no existing distRun
     AND ((chipBackgroundRun.state = 'full') OR (distTarget.clean AND chipBackgroundRun.state = 'cleaned'))

branches/eam_branches/ipp-20140402/ippTools/share/disttool_definebyquery_diff.sql

@@ -25 +25 @@
 WHERE  distTarget.state = 'enabled'
     AND rcInterest.state = 'enabled'
-    AND distRun.dist_id IS NULL
     AND ((diffRun.state = 'full') OR (distTarget.clean AND diffRun.state = 'cleaned'))
 

branches/eam_branches/ipp-20140402/ippTools/share/disttool_definebyquery_fake.sql

@@ -20 +20 @@
 WHERE  distTarget.state = 'enabled'
     AND rcInterest.state = 'enabled'
-    AND distRun.dist_id IS NULL
     AND ((fakeRun.state = 'full') OR (distTarget.clean AND fakeRun.state = 'cleaned'))

branches/eam_branches/ipp-20140402/ippTools/share/disttool_definebyquery_ff.sql

@@ -21 +21 @@
 WHERE  distTarget.state = 'enabled'
     AND rcInterest.state = 'enabled'
-    AND distRun.dist_id IS NULL
     AND fullForceRun.state = 'full'
     AND stackRun.filter = distTarget.filter

branches/eam_branches/ipp-20140402/ippTools/share/disttool_definebyquery_raw.sql

@@ -16 +16 @@
 JOIN rcInterest USING(target_id)
 LEFT JOIN distRun ON distRun.stage_id = exp_id AND distRun.target_id = distTarget.target_id
-WHERE distRun.dist_id IS NULL           -- no existing distRun for this exposure
-    AND distTarget.state = 'enabled'    -- target and intrest are enabled
+WHERE distTarget.state = 'enabled'    -- target and intrest are enabled
     AND rcInterest.state = 'enabled'
     AND magicDSRun.state = 'full'       -- destreaked files are available

branches/eam_branches/ipp-20140402/ippTools/share/disttool_definebyquery_sky.sql

@@ -23 +23 @@
     AND rcInterest.state = 'enabled'
     AND distTarget.filter = 'multi'
-    AND distRun.dist_id IS NULL
     AND ((staticskyRun.state = 'full') OR (distTarget.clean AND staticskyRun.state = 'cleaned'))
     -- we shouldn't need to check fault. If faulted it shouldn't be full

branches/eam_branches/ipp-20140402/ippTools/share/disttool_definebyquery_sky_singlefilter.sql

@@ -22 +22 @@
 WHERE  distTarget.state = 'enabled'
     AND rcInterest.state = 'enabled'
-    AND distRun.dist_id IS NULL
     AND ((staticskyRun.state = 'full') OR (distTarget.clean AND staticskyRun.state = 'cleaned'))
     AND stackRun.filter = distTarget.filter

branches/eam_branches/ipp-20140402/ippTools/share/disttool_definebyquery_skycal.sql

@@ -1 @@
-SELECT
+SELECT DISTINCT
     'skycal' as stage,
     skycalRun.skycal_id AS stage_id,
@@ -21 +21 @@
 WHERE  distTarget.state = 'enabled'
     AND rcInterest.state = 'enabled'
-    AND distRun.dist_id IS NULL
     AND ((skycalRun.state = 'full') OR (distTarget.clean AND skycalRun.state = 'cleaned'))
     AND stackRun.filter = distTarget.filter

branches/eam_branches/ipp-20140402/ippTools/share/disttool_definebyquery_stack.sql

@@ -21 +21 @@
 WHERE  distTarget.state = 'enabled'
     AND rcInterest.state = 'enabled'
-    AND distRun.dist_id IS NULL
     AND ((stackRun.state = 'full') OR (distTarget.clean AND stackRun.state = 'cleaned'))
     -- we shouldn't need to check fault. If faulted it shouldn't be full

branches/eam_branches/ipp-20140402/ippTools/share/disttool_definebyquery_warp.sql

@@ -22 +22 @@
 WHERE  distTarget.state = 'enabled'
     AND rcInterest.state = 'enabled'
-    AND distRun.dist_id IS NULL
     AND ((warpRun.state = 'full') OR (distTarget.clean AND warpRun.state = 'cleaned'))
 

branches/eam_branches/ipp-20140402/ippTools/share/disttool_definebyquery_warp_bg.sql

@@ -24 +24 @@
 WHERE  distTarget.state = 'enabled'
     AND rcInterest.state = 'enabled'
-    AND distRun.dist_id IS NULL
     AND ((warpBackgroundRun.state = 'full') OR (distTarget.clean AND warpBackgroundRun.state = 'cleaned'))
 

branches/eam_branches/ipp-20140402/ippTools/share/warptool_warped.sql

@@ -4 +4 @@
     warpRun.workdir,
     warpRun.label,
+    warpRun.data_group,
     warpImfile.warp_skyfile_id,
     rawExp.exp_id,

branches/eam_branches/ipp-20140402/ippTools/src/fftool.c

@@ -498 +498 @@
     PXOPT_COPY_STR(config->args, where, "-filter",      "stackRun.filter", "LIKE");
     PXOPT_COPY_S16(config->args, where, "-fault",      "fullForceResult.fault", "==");
+    PXOPT_COPY_S16(config->args, where, "-quality",    "fullForceResult.quality", "==");
     pxskycellAddWhere(config, where);
 

branches/eam_branches/ipp-20140402/ippTools/src/fftoolConfig.c

@@ -116 +116 @@
     psMetadataAddStr(resultArgs, PS_LIST_TAIL, "-label", 0, "search by label", NULL);
     psMetadataAddStr(resultArgs, PS_LIST_TAIL, "-data_group", 0, "search by data_group (LIKE comparison)", NULL);
+    psMetadataAddS16(resultArgs, PS_LIST_TAIL, "-quality", 0, "search by quality value", 0);
     pxskycellAddArguments(resultArgs);
     psMetadataAddS16(resultArgs, PS_LIST_TAIL, "-fault", 0, "search by fault code", 0);

branches/eam_branches/ipp-20140402/ippconfig/gpc1/ppStack.config

@@ -167 +167 @@
     PSF.TARGET.AS.MAX  BOOL  TRUE
     PSF.TARGET.AS.MAX.EPSILON F32 0.1
+    THRESHOLD.MASK  F32     0.01             # Threshold for mask deconvolution (0..1)
+    COMBINE.ITER    F32     0.5             # Number of rejection iterations per input
+    COMBINE.REJ     F32     2.5             # Rejection threshold in combination (sigma)
+    COMBINE.SYS     F32     0.1             # Relative systematic error in combination
+    COMBINE.DISCARD F32     0.2             # Discard fraction for Olympic weighted mean
+    IMAGE.REJ       F32     0.1             # Rejected pixel fraction threshold for rejecting entire image
 END
 
@@ -183 +189 @@
     PSF.TARGET.AS.MAX  BOOL  TRUE
     PSF.TARGET.AS.MAX.EPSILON F32 0.1
+    THRESHOLD.MASK  F32     0.5             # Threshold for mask deconvolution (0..1)
+    COMBINE.ITER    F32     0.5             # Number of rejection iterations per input
+    COMBINE.REJ     F32     3.5             # Rejection threshold in combination (sigma)
+    COMBINE.SYS     F32     0.1             # Relative systematic error in combination
+    COMBINE.DISCARD F32     0.2             # Discard fraction for Olympic weighted mean
+    IMAGE.REJ       F32     0.1             # Rejected pixel fraction threshold for rejecting entire image
+END
+
+STACK_MDTHREEPI_DEEP_1DG     METADATA   ## revised config for MD+3PI 
+    OUTPUT.NOCOMP           BOOL  TRUE
+    OUTPUT.LOGFLUX          BOOL  FALSE
+    OUTPUT.REPLICATE        BOOL  TRUE
+    STACK.TYPE              STR   DEEP_STACK
+    PSF.INPUT.CLIP.SIMPLE   BOOL  True
+    PSF.INPUT.MAX           F32   10.0    #
+    PSF.INPUT.CLIP.NSIGMA   F32   100.0
+    PSF.INPUT.THRESH        F32   NAN
+    PSF.INPUT.ASYMMETRY     F32   NAN
+    MATCH.REJ               F32   4.0     # need ot tune?
+    SAFE                    BOOL  FALSE
+    PSF.TARGET.AS.MAX  BOOL  TRUE
+    PSF.TARGET.AS.MAX.EPSILON F32 0.1
+    THRESHOLD.MASK  F32     0.01             # Threshold for mask deconvolution (0..1)
+    COMBINE.ITER    F32     0.5             # Number of rejection iterations per input
+    COMBINE.REJ     F32     2.5             # Rejection threshold in combination (sigma)
+    COMBINE.SYS     F32     0.1             # Relative systematic error in combination
+    COMBINE.DISCARD F32     0.2             # Discard fraction for Olympic weighted mean
+    IMAGE.REJ       F32     0.1             # Rejected pixel fraction threshold for rejecting entire image
 END
 

branches/eam_branches/ipp-20140402/ippconfig/recipes/nightly_science.config

@@ -3 +3 @@
   MODE           STR CHIP
   COMMAND        STR chiptool -dbname @DBNAME@ -updaterun -set_state goto_cleaned -state full -set_label goto_cleaned -label @LABEL@ -data_group @DATA_GROUP@
-  RETENTION_TIME U16 7
+  RETENTION_TIME U16 5 
 END
 CLEAN_MODES METADATA
@@ -9 +9 @@
   COMMAND        STR warptool -dbname @DBNAME@ -updaterun -set_state goto_cleaned -state full -set_label goto_cleaned -label @LABEL@ -data_group @DATA_GROUP@
 #  COMMAND        STR true
-  RETENTION_TIME U16 7
+  RETENTION_TIME U16 5
 END
 CLEAN_MODES METADATA
   MODE           STR DIFF
   COMMAND        STR difftool -dbname @DBNAME@ -updaterun -set_state goto_cleaned -state full -set_label goto_cleaned -label @LABEL@ -data_group @DATA_GROUP@
-  RETENTION_TIME S16 7
+  RETENTION_TIME S16 5
 END
 CLEAN_MODES METADATA

branches/eam_branches/ipp-20140402/ippconfig/recipes/ppStack.config

@@ -189 +189 @@
 END
 
+STACK_MDTHREEPI_DEEP_1DG   METADATA
+END
+
 
 

branches/eam_branches/ipp-20140402/ippconfig/recipes/psphot.config

@@ -491 +491 @@
   PEAKS_NMAX_TOTAL                    S32   0 # set this to limit the allowed number of peaks - Yields fault instead of avoid memory explosion
 
-  EXT_MODEL_TYPE_FOR_CFF              STR   PS_MODEL_DEV
+#  EXT_MODEL_TYPE_FOR_CFF              STR   PS_MODEL_DEV
 
   SAVE.RESID                          BOOL  TRUE

branches/eam_branches/ipp-20140402/ippconfig/recipes/reductionClasses.mdc

@@ -663 +663 @@
       STACKPHOT_SINGLE_PSPHOT STR     STACKPHOT_SINGLE
       STATICSKY_CALIBRATION   STR     STATICSKY_DEEPCAL 
+END
+
+MDTHREEPI_DEEP_STACK_1DG            METADATA
+      STACK_PPSTACK           STR     STACK_MDTHREEPI_DEEP_1DG
+      STACK_PPSUB             STR     STACK_1DG
+      STACK_PSPHOT            STR     STACK
+      STACKPHOT_PSPHOT        STR     STACKPHOT
+      STACKPHOT_PPSUB         STR     STACKPHOT_1DG
+      STACKPHOT_PPSTACK       STR     STACKPHOT
+      STACKPHOT_SINGLE_PSPHOT STR     STACKPHOT_SINGLE
+      STATICSKY_CALIBRATION   STR     STATICSKY_DEEPCAL
 END
 

branches/eam_branches/ipp-20140402/ppBackground/src/ppBackgroundStack.h

@@ -23 +23 @@
 
   psImageMap *modelMap;
+  psS32 model_iteration;
   // These are the full extent of the input data
   psF32 ra_min;

branches/eam_branches/ipp-20140402/ppBackground/src/ppBackgroundStackCamera.c

@@ -114 +114 @@
         }
       }
-
-/*       psMemStats(0,&A,&P); */
-/*       fprintf(stderr,"fpa %ld %ld\n",A,P); */
+      
+      psF32 exptime = 1.0;
+      exptime = psMetadataLookupF32(NULL, phu->header, "EXPTIME");
 
       pmChip *chip;                       // Chip from FPA
@@ -197 +197 @@
             decim->data.F32[v][u] = tp->y;
             model->data.F32[v][u] = 0.0;
-
+            image->data.F32[v][u] /= exptime;
             // Check the bounds so we'll know how large of an area to model in the map
             if (tp->x < data->ra_min) { data->ra_min = tp->x; }
@@ -277 +277 @@
           }
         }
-
-        
-/*      psMemStats(0,&A,&P); */
-/*      fprintf(stderr,"chip2 %ld %ld\n",A,P); */
       } // end chip loop
 
@@ -290 +286 @@
                        "N", PS_META_REPLACE,
                        "counter", i);
-      // data->models->name->XY__->{image/ra/dec/calibrated/offset/scale}
-      // data->models->counter->XY__->{image/ra/dec/calibrated/offset/scale}
-
-      
       
       if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {

branches/eam_branches/ipp-20140402/ppBackground/src/ppBackgroundStackData.c

@@ -50 +50 @@
 
     data->modelMap = NULL;
+    data->model_iteration = 0;
     data->ra_min   = 1e9;
     data->ra_max   = -1e9;

branches/eam_branches/ipp-20140402/ppBackground/src/ppBackgroundStackLoop.c

@@ -9 +9 @@
 
 #include "ppBackgroundStack.h"
+
+#define WCS_TOLERANCE 0.001             // Tolerance for WCS
+
 
 bool ppBackgroundStackLoop(ppBackgroundStackData *data // Run-time data
@@ -76 +79 @@
       pmFPAview *view = pmFPAviewAlloc(0);
 
+      //      pmHDU *phu = pmFPAviewThisPHU(view, stack->fpa);
+      psF32 exptime = 1.0;
+      exptime = psMetadataLookupF32(NULL, stack->fpa->hdu->header, "EXPTIME");
+      
       // PART 1:
       // Determine the extent of the model map for this stack
@@ -116 +123 @@
       if (tp->y > data->y_max) { data->y_max = tp->y; }
 
-/*       data->x_min -= data->ra_min; */
-/*       data->x_max -= data->ra_min; */
-/*       data->y_min -= data->dec_min; */
-/*       data->y_max -= data->dec_min; */
-      
-      
       psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
       psImageBinning *binning = psImageBinningAlloc();
-      binning->nXruff = 15; // Number of samples
-      binning->nYruff = 15; 
-      //    binning->nXfine = ceil(data->ra_max - data->ra_min) + 1; // This is the range we're looking at
-      //    binning->nYfine = ceil(data->dec_max - data->dec_min) + 1;
+      binning->nXruff = 13; // Number of samples
+      binning->nYruff = 13; 
       binning->nXfine = ceil(data->x_max - data->x_min) + 1;
       binning->nYfine = ceil(data->y_max - data->y_min) + 1;
@@ -142 +141 @@
       P_PSIMAGE_SET_ROW0(sizeImage, data->y_min);
       data->modelMap = psImageMapAlloc(sizeImage,binning,stats);
-/*       psFree(sizeImage); */
-/*       data->modelMap = psImageMapNoImageAlloc( binning,stats); */
-/*       P_PSIMAGE_SET_COL0(data->modelMap->map, (data->x_min - binning->nXskip) / binning->nXbin); */
-/*       P_PSIMAGE_SET_ROW0(data->modelMap->map, (data->y_min - binning->nYskip) / binning->nYbin); */
-
-      // force col0/row0
-/*       data->modelMap->map->col0 = data->modelMap->binning->nXskip; */
-/*       data->modelMap->map->row0 = data->modelMap->binning->nYskip; */
       
       // PART 2:
@@ -169 +160 @@
 
       // This is where an iterative solution loop would likely start.
-      {
+      for (int iterator = 0; iterator < 4; iterator++) {
         // Construct the offset information
         printf("Model fit!\n");
@@ -298 +289 @@
                 psPlaneTransformApply(tp, stack->fpa->toTPA, pix);
                 
-                model->image->data.F32[y][x] = psImageMapEval(data->modelMap,tp->x,tp->y);
+                model->image->data.F32[y][x] = exptime * psImageMapEval(data->modelMap,tp->x,tp->y);
                 corr->image->data.F32[y][x] = readout->image->data.F32[y][x] + model->image->data.F32[y][x];
                 
@@ -306 +297 @@
             psFree(fp);
             psFree(tp);
+
+            // Copy WCS (from ppStackUpdateHeader)
+            pmHDU *inHDU = pmHDUFromCell(readout->parent);
+            model->parent->hdu = pmHDUAlloc(NULL);
+            corr->parent->hdu = pmHDUAlloc(NULL);
+            pmHDU *modHDU= pmHDUFromCell(model->parent);
+            pmHDU *corHDU= pmHDUFromCell(corr->parent);
+
+            if (!modHDU || !inHDU) {
+              psWarning("Unable to find HDU at FPA level to copy wcs!");
+            }
+            else {
+              if (!pmAstromReadWCS(stack_model->fpa,model_cell->parent,inHDU->header,1.0)) {
+                psErrorClear();
+                psWarning("Unable to read WCS astrometry from input FPA!");
+              }
+              else {
+                if (!modHDU->header) {
+                  modHDU->header = psMetadataAlloc();
+                }
+                if (!pmAstromWriteWCS(modHDU->header, stack_model->fpa,model_cell->parent, WCS_TOLERANCE)) {
+                  psErrorClear();
+                  psWarning("Unable to read WCS astrometry from input FPA!");
+                }
+              }
+              if (!pmAstromReadWCS(stack_corr->fpa,corr_cell->parent,inHDU->header,1.0)) {
+                psErrorClear();
+                psWarning("Unable to read WCS astrometry from input FPA!");
+              }
+              else {
+                if (!corHDU->header) {
+                  corHDU->header = psMetadataAlloc();
+                }
+                if (!pmAstromWriteWCS(corHDU->header, stack_corr->fpa,corr_cell->parent, WCS_TOLERANCE)) {
+                  psErrorClear();
+                  psWarning("Unable to read WCS astrometry from input FPA!");
+                }
+              }
+            } // End WCS saving.
+
+            
           } // Close readout
           printf("    I'm done with that readout\n");
@@ -335 +367 @@
       psFree(view);
       psFree(data->modelMap);
+      psFree(sizeImage);
     }
                 

branches/eam_branches/ipp-20140402/ppBackground/src/ppBackgroundStackMath.c

@@ -55 +55 @@
           }
           psMetadata *chipData = psMetadataLookupPtr(NULL, expItem->data.md, workingChip);
+          if (!chipData) { continue; }
           psImage *image      = psMetadataLookupPtr(NULL, chipData, "bkg image");
+          if (!image) { continue; }
           psVectorAppend(tmp,image->data.F32[v][u]);
         } // End loop over exposures
@@ -63 +65 @@
 
         psFree(expIter);
+        psFree(tmp);
       } // End u
     } // End v
+    // Remove the median value from this data.  We just want the tilts, not the offsets
+    psStatsInit(stats);
+    psImageStats(stats,solution,NULL,0);
+    for (v = 0; v < solution->numRows; v++) {
+      for (u = 0; u < solution->numCols; u++) {
+        solution->data.F32[v][u] -= stats->robustMedian;
+      }
+    }
+    
+    
   } // End working chip scan
 
@@ -90 +103 @@
     psMetadataItem *chipItem;
     while ((chipItem = psMetadataGetAndIncrement(chipIter))) {
+      //      const char *chipName = chipItem->name;
+      
       psImage *image = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg image");
       psImage *ra    = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra");
       psImage *dec   = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec");
-      
+      //      psImage *camera= psMetadataLookupPtr(NULL, data->OTA_solutions, chipName);      
       psVector *obs  = psVectorAllocEmpty(image->numRows, PS_TYPE_F32);
       psVector *model= psVectorAllocEmpty(image->numRows, PS_TYPE_F32);
-      
+
+      int j = 0;
+      int used = 0;
       for (v = 0; v < image->numRows; v++) {
         for (u = 0; u < image->numCols; u++) {
           if ((ra->data.F32[v][u] < data->x_min)||(ra->data.F32[v][u] > data->x_max)||
-              (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) { continue; }
-          psVectorAppend(obs,image->data.F32[v][u]);
+              (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) { j++; continue; }
+          psVectorAppend(obs,image->data.F32[v][u]);// - camera->data.F32[v][u]);
           psVectorAppend(model, psImageMapEval(data->modelMap,ra->data.F32[v][u],dec->data.F32[v][u]));
-          //      obs->data.F32[v * image->numCols + u] = image->data.F32[v][u];
-          //      model->data.F32[v * image->numCols + u] = psImageMapEval(data->modelMap,ra->data.F32[v][u],dec->data.F32[v][u]);
-        }
-      }
-
-      psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD,1);
-      int status = psVectorFitPolynomial1D(poly,NULL,0,model,NULL,obs);
-      if (!status) {
-        psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg offset", PS_META_REPLACE, "background offset for this exposure/ota pair", poly->coeff[0]);
-        psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg scale", PS_META_REPLACE, "background scale for this exposure/ota pair", poly->coeff[1]);
-      }
-      psFree(poly);
+          j++;
+          used++;
+        }
+      }
+
+      if (used > 0) {
+
+        psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD,1);
+        int status = psVectorFitPolynomial1D(poly,NULL,0,model,NULL,obs);
+        printf("in model fit loop: %d %d %d %f %f\n",status,j,used,poly->coeff[0],poly->coeff[1]);
+        if (status && (poly->coeff[1] != 0.0)) {
+          psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg offset", PS_META_REPLACE, "background offset for this exposure/ota pair", poly->coeff[0]);
+          psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg scale", PS_META_REPLACE, "background scale for this exposure/ota pair", poly->coeff[1]);
+        }
+        psFree(poly);
+      }
     } // End OTA loop
     psFree(chipIter);
@@ -155 +176 @@
               (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) { continue; }
 
-          model->data.F32[v][u] = scale * image->data.F32[v][u] - offset - camera->data.F32[v][u];
+          model->data.F32[v][u] = scale * image->data.F32[v][u] + offset - camera->data.F32[v][u];
         }
       }
@@ -170 +191 @@
 // This "averaging" is done using the psImageMapClipFit.
 bool ppBackgroundStackModelFit(ppBackgroundStackData *data) {
-  long j;
-  int u,v;
-
-  long used;
+  long j = 0;
+  int u,v;
+
+  long used = 0;
   psS16 N = psMetadataLookupS16(NULL, data->models, "N");
   psVector *X = psVectorAllocEmpty(N * 13 * 13,PS_TYPE_F32);
@@ -197 +218 @@
       psImage *ra    = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra");
       psImage *dec   = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec");
-      //      psImage *model = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg image");
+#define DUMP_DATA 0
+#if DUMP_DATA
+      psImage *model = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg image");
+
+      psF32 offset = psMetadataLookupF32(NULL,chipItem->data.md,"bkg offset");
+      psF32 scale  = psMetadataLookupF32(NULL,chipItem->data.md,"bkg scale");
+
+#endif
       for (v = 0; v < calib->numRows; v++) {
         for (u = 0; u < calib->numCols; u++) {
@@ -216 +244 @@
           used++;
           j++;
-          
-/*        printf("DATA %ld %ld %f %f %f %f\n", */
-/*               j,used, */
-/*               ra->data.F32[v][u], */
-/*               dec->data.F32[v][u], */
-/*               calib->data.F32[v][u], */
-/*               model->data.F32[v][u]); */
-                 
+#if DUMP_DATA
+          printf("DATA %d %ld %ld %f %f %f %f %f %f\n",
+                 data->model_iteration,j,used,
+                 offset,scale,
+                 ra->data.F32[v][u],
+                 dec->data.F32[v][u],
+                 calib->data.F32[v][u],
+                 model->data.F32[v][u]);
+#endif 
         }
       }
@@ -232 +261 @@
   bool fitStatus;
   bool status = psImageMapClipFit(&fitStatus,data->modelMap,stats, mask, 1, X, Y, Z, E);
-
+  data->model_iteration++;
   psFree(expIter);
   psFree(X);

branches/eam_branches/ipp-20140402/ppSkycell/src/ppSkycellLoop.c

@@ -125 +125 @@
         for (int j = 0; j < numProj && !found; j++) {
             if (wcs->crval1 == crval1->data.F64[j] && wcs->crval2 == crval2->data.F64[j] &&
-                wcs->cdelt1 == cdelt1->data.F64[j] && wcs->cdelt2 == cdelt1->data.F64[j]) {
+                wcs->cdelt1 == cdelt1->data.F64[j] && wcs->cdelt2 == cdelt2->data.F64[j]) {
               //                regionMinMax(projRegions->data[j], region);
               psRegion *proj = projRegions->data[j];

branches/eam_branches/ipp-20140402/psModules/src/imcombine/pmStack.c

@@ -43 +43 @@
 // #define TEST_X 972
 // #define TEST_Y 3213
-#define TEST_X 3289
-#define TEST_Y 4810
-#define TEST_RADIUS 0.5                 // Radius to examine
+//#define TEST_X 3289
+//#define TEST_Y 4810
+//#define TEST_RADIUS 0.5                 // Radius to examine
+//MEH -- streak-like junk md04s065i
+#define TEST_X 1129 
+#define TEST_Y 4256
+#define TEST_RADIUS 2.0                 // Radius to examine
 # endif
 
@@ -108 +112 @@
 
 // KMM functions to do bimodality rejection of pixels
-float gaussian(float x, float m, float s) {
+double gaussian(float x, float m, float s) {
   return(pow(s * sqrt(2 * M_PI),-1) * exp(-0.5 * pow( (x - m) / s, 2)));
 }
@@ -116 +120 @@
                          float *mU, float *sU,
                          float *pi1, float *m1, float *s1,
-                         float *pi2, float *m2, float *s2) {
+                         float *pi2, float *m2, float *s2,
+                         int xyrdebug) {
   assert(values);
   assert(values->type.type == PS_TYPE_F32);
@@ -151 +156 @@
   }
 
+  if (xyrdebug == 1) { 
+      fprintf(stderr,"KMM uni: %d %f %d (%f %f)\n", 
+              xyrdebug,logL_unimodal,discrepant_index, 
+              *mU,*sU); 
+  } 
+
   // Do EM loop
   float dL = 0;
@@ -186 +197 @@
   *pi2 = 0.5;
 
-  float g1,g2,norm;
+  //MEH -- need to be double to help avoid 0 in norm
+  double g1,g2,norm;
   float w1,w2;
 
@@ -203 +215 @@
 /*            *m2,*s2,*pi2); */
 /*     } */
+
+    if (xyrdebug == 1) { 
+        fprintf(stderr,"KMM EM iter: %d %f %f %f %f (%f %f %e) (%f %f %e)\n", 
+                *iter,logL_unimodal,logL_bimodal,oldL,dL, 
+                *m1,*s1,*pi1, 
+                *m2,*s2,*pi2); 
+    } 
+
     // Expectation/P-stage
     for (i = 0; i < values->n; i++) { // Calculate probabilities for each mode
@@ -208 +228 @@
       g2 = gaussian(values->data.F32[i],*m2,*s2);
       norm = (*pi1 * g1 + *pi2 * g2);
-      P1->data.F32[i] = (*pi1 * g1) / norm;
-      P2->data.F32[i] = (*pi2 * g2) / norm;
+      //MEH -- must protect denom from norm=0
+      if (norm > 0) {
+          P1->data.F32[i] = (*pi1 * g1) / norm;
+          P2->data.F32[i] = (*pi2 * g2) / norm;
+      } else {
+          P1->data.F32[i] = 0.0;
+          P2->data.F32[i] = 0.0;
+      }  
+ 
+      if (xyrdebug == 1) {
+          fprintf(stderr,"KMM EM-P loop: %d %d %le %le %le\n",
+                         *iter,i,norm,g1,g2);
+      }
+
     }
     // Maximization/M-stage
@@ -231 +263 @@
       w1 += P1->data.F32[i];
       w2 += P2->data.F32[i];
+
+      if (xyrdebug == 1) {
+          fprintf(stderr,"KMM EM-M loop: %d %d (%f %f %f %e) (%f %f %f %e)\n",
+                         *iter,i,*m1,values->data.F32[i],w1,P1->data.F32[i],*m2,values->data.F32[i],w2,P2->data.F32[i]);
+      }
+
     }
     *m1 /= w1;
@@ -250 +288 @@
       *pi2 = 0.0;
     }
-    if (*s1 == 0) { // sigma may not be zero
-      *s1 = KMM_SMALL_NUMBER * *m1;
-    }
-    if (*s2 == 0) { // sigma may not be zero
-      *s2 = KMM_SMALL_NUMBER * *m2;
-    }
+    if (*s1 == 0) { // sigma may not be zero -- MEH -- nor <0 and need additive offset if m~0
+      *s1 = fabsf(KMM_SMALL_NUMBER * *m1) + KMM_SMALL_NUMBER;
+    }
+    if (*s2 == 0) { // sigma may not be zero 
+      *s2 = fabsf(KMM_SMALL_NUMBER * *m2) + KMM_SMALL_NUMBER;
+    }
+
+    if (xyrdebug == 1) { 
+        fprintf(stderr,"KMM EM end: %d %f %f %f %f (%f %e %e %f) (%f %e %e %f)\n", 
+                *iter,logL_unimodal,logL_bimodal,oldL,dL, 
+                *m1,*s1,*pi1,w1, 
+                *m2,*s2,*pi2,w2); 
+    } 
+
   } // End EM phase
 
@@ -267 +313 @@
     *Punimodal = 1.0;
   }
+
+  if (xyrdebug == 1) { 
+      fprintf(stderr,"KMM calc Puni: %d %f %d %f\n", 
+              xyrdebug,lambda,df,*Punimodal); 
+  } 
+
   psFree(P1);
   psFree(P2);
 }
 
-static void KMMFindPopular(const psVector *values, float *Punimodal, float *mean, float *sigma, float *pi) {
+static void KMMFindPopular(const psVector *values, float *Punimodal, float *mean, float *sigma, float *pi, int xyrdebug) {
   float KMM_MINIMUM_PVALUE = 0.05; // Should be an option.
   float mU,sU;
@@ -283 +335 @@
                &mU,&sU,
                &pi1,&m1,&s1,
-               &pi2,&m2,&s2);
+               &pi2,&m2,&s2,xyrdebug);
 /*   fprintf(stdout,"%g %g : %g %g %g : %g %g %g : %g %d\t", */
 /*        mU,sU, */
@@ -898 +950 @@
     psVector *pixelSuspects = buffer->suspects; // Is the pixel suspect?
     psVector *pixelLimits = buffer->limits; // Is the pixel suspect?
+    //MEH -- adding a debug option for TESTING xyr position but could be better..
+    int xyrdebug = 0;
 
     // KMM values;
@@ -917 +971 @@
         // This should probably be an option
         if (useKMM) {
-          KMMFindPopular(pixelData,&Punimodal,&KMMmean,&KMMsigma,&KMMpi);
+#ifdef TESTING
+            if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
+                xyrdebug = 1;
+            }
+# endif
+          KMMFindPopular(pixelData,&Punimodal,&KMMmean,&KMMsigma,&KMMpi,xyrdebug);
           CHECKPIX(x,y,"KMM Popularity Contest: (%d,%d) Puni: %g Mean: %f Sigma %f Pi: %f\n",
                    x,y,Punimodal,KMMmean,KMMsigma,KMMpi);

branches/eam_branches/ipp-20140402/psModules/src/imcombine/pmSubtraction.c

@@ -37 +37 @@
 #define USE_KERNEL_ERR                  // Use kernel error image?
 #define NUM_COVAR_POS 5                 // Number of positions for covariance calculation
+//MEH -- this is causing diffim fault 5 -- seems not as robust
 #define USE_LOGFIT_REJECT
 

branches/eam_branches/ipp-20140402/psModules/src/objects/pmPSFtryFitPSF.c

@@ -125 +125 @@
 
     if (Npsf == 0) {
+#if 0
         // DEBUG code: save the PSF model fit data in detail
 
@@ -143 +144 @@
 
             pmSource *source = psfTry->sources->data[i];
+
+            if (!source->modelPSF) continue;
 
             float par7 = (source->modelPSF->params->n == 7) ? -100 : source->modelPSF->params->data.F32[PM_PAR_7];
@@ -155 +158 @@
         }
         fclose (f);
-
+#endif
         psError(PS_ERR_UNKNOWN, false, "No sources with good PSF fits after model is built.");
         return false;

branches/eam_branches/ipp-20140402/psastro/src/psastroLoadCrosstalk.c

@@ -108 +108 @@
   SPIKE_MAX_MAG = psMetadataLookupF32(&status, recipe, "REFSTAR_MASK_BLEED_MAG_MAX");
   SPIKE_MAX_MAG += MagOffset;
-  
+  //  psTraceSetLevel("psastro.crosstalk",2);  
   psTrace("psastro.crosstalk",2,"%f %f %f %f\n",zeropt,exptime,MAX_MAG,MagOffset);
 
@@ -147 +147 @@
           // Check each reference star
           for (int i = 0; i < refstars->n; i++) {
+            // This is the source of the ct.
             pmAstromObj *ref = refstars->data[i];
 
@@ -174 +175 @@
 
             float x_t_chip,y_t_chip;
-
+            int faint_ct = 0;
             pmCellCoordsForChip(&x_cell,&y_cell,CTsourceCell,ref->chip->x,ref->chip->y);
 
@@ -180 +181 @@
             y_t_cell = y_cell;
 
+            // 2014-04-15 CZW: Here are the list of known cross talks:
+            // Inter-chip
+            // OTA2yXY3v <=> OTA3yXY3v
+            // OTA4yXY3v <=> OTA5yXY3v
+            // Intra-chip
+            // OTA2yXY5v <=> OTA2yXY6v
+            // OTA5yXY5v <=> OTA5yXY6v
+            // One way fainter
+            // OTA2yXY7v => OTA3yXY2v
+            // OTA5yXY7v => OTA4yXY2v
             // Determine if this combination of chip and cell produces a crosstalk artifact
-            if ((X == 2)&&(! ((U == 3)||(U == 5)||(U == 6)))) {
+            if ((X == 2)&&(! ((U == 3)||(U == 5)||(U == 6)||(U == 7)))) {
               psTrace ("psastro.crosstalk",2,"Cell (%d%d) on chip (%d%d) not a known crosstalk source.",U,V,X,Y);
+              psTrace ("psastro.crosstalk",2,"t1: %d t2: %d",
+                       (X == 2),
+                       (! ((U == 3)||(U == 5)||(U == 6)||(U == 7))));
               continue;
             }
@@ -193 +207 @@
               continue;
             }
-            if ((X == 5)&&(! ((U == 3)||(U == 5)||(U == 6)))) {
+            if ((X == 5)&&(! ((U == 3)||(U == 5)||(U == 6)||(U == 7)))) {
               psTrace ("psastro.crosstalk",2,"Cell (%d%d) on chip (%d%d) not a known crosstalk source.",U,V,X,Y);
               continue;
@@ -212 +226 @@
                 Ut = 5;
               }
+              if (U == 7) {
+                Xt = 3;
+                Ut = 2;
+                faint_ct = 1;
+              }
             }
             else if (X == 3) {
@@ -229 +248 @@
                 Ut = 5;
               }
+              if (U == 7) {
+                Xt = 4;
+                Ut = 2;
+                faint_ct = 1;
+              }
             }
 
@@ -284 +308 @@
 
             // Determine if we need to add a spike mask.
+            // CZW 2014-04-15 I'm not adding the new faint CT to this operation.  There's no evidence they produce such things.
 /*          psWarning("BLAH: %g %g %g\n",ref->Mag,MagOffset,SPIKE_MAX_MAG); */
             if (ref->Mag < SPIKE_MAX_MAG) {

branches/eam_branches/ipp-20140402/psphot/src/psphotFullForceReadout.c

@@ -78 +78 @@
     psphotMergeSources (config, view, filerule);
 
-
     // generate a psf model for any readouts which need one
     // psphotFullForcePSF (config, view, filerule);
@@ -84 +83 @@
     // sources; we should mark with a flag bit the ones we actually want to use as PSF
     // stars (this means we need to supply this info in the load).
-    psphotChoosePSF (config, view, filerule, false);
+    if (!psphotChoosePSF (config, view, filerule, false)) {
+        // PSPHOT_ERR_DATA causes this program to exit gracefully
+        psError (PSPHOT_ERR_DATA, false, "failed to construct psf model");
+        return  psphotReadoutCleanup (config, view, filerule);
+    }
 
     // Construct an initial model for each object, set the radius to fitRadius, set circular

branches/eam_branches/ipp-20140402/psphot/src/psphotFullForceSummary.c

@@ -193 +193 @@
         if (ID >= sortedTable->n) {
             sortedTable = psArrayRealloc(sortedTable, 2*ID);
+            // Why doesn't psArrayRealloc do this?????
+            sortedTable->n = sortedTable->nalloc;
         }
         if (sortedTable->data[ID]) {

branches/eam_branches/ipp-20140402/psphot/src/psphotLensing.c

@@ -193 +193 @@
     // set the max order (0 = constant) which the number of psf stars can support:
     // we require only 3 stars for n = 0, increase stars / cell for higher order
+    if (Npsf < 3) {
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources passed selection criteria, skipping lensing psf trends");
+        return true;
+    }
     int LensingTrendOrder = 0;
     if (Npsf >=  16) LensingTrendOrder = 1; // 4 cells

branches/eam_branches/ipp-20140402/tools/dvodb_calib/README

@@ -9 +9 @@
 Step Four: run gene.relastro.part2
 
-Step Five: run gene.relphot
+Step Five: run gene.relphot.3pinew or run gene.relphot.mdnew (depends on tess)