/***************************************************************************** ****************************************************************************** ** ** ** FILE: tclMerge.cc ** ** ** This file contains TCL functions to merge, unmerge, and generate ** merging statistics. ** ** ** ** ** ENVIRONMENT: ** C++. ** ** REQUIRED PRODUCTS: ** photoSch ** ****************************************************************************** ****************************************************************************** */ #include #include #include #include #include #include #include #include #include #include #include "tsMerge.h" #include #include extern "C" { #include #include } #include "arObjectStatus.h" #include #include #include #include #include #include #include #include const int nFilters = 5; /* Number of filters */ const int nColors = 4; /* Number of colors */ /*============================================================================ ** ** ** TCL VERB: xSetObjectStatusDontCall ** ** ** Set the object status for all objects in a single Frames Pipeline Run, ** and merge objects in the overlap regions between adjacent fields. ** This consists of three steps, performed in the following order: **
    **
  1. Set each object's status (i.e., set/don't set the GOOD flag) ** based on the values of the object flags. **
  2. Merge objects lying in the overlap regions between adjacent ** fields. **
  3. Set/don't set the OK_RUN and DUPLICATE (for GOOD objects only), ** determining whether this object is a primary detection for this run or ** not. This includes resolving the matched pairs of objects in the ** overlap region, requiring that exactly one object in each matched pair ** has its status set to OK_RUN. **
** The following object flags are set by this routine: **
**
SET **
All objects will have this flag set, indicating the objects have ** been processed by "setObjectStatus". **
GOOD **
If set, then this object has been properly measured, as ** determined by the object flags, to the extent ** that it MAY be used in the primary survey, depending on the status of ** its other flags. If not set, then the object has NOT been properly ** measured, and it may NOT be used in the primary survey. **
** For each GOOD object, the following two flags MAY be set: **
**
OK_RUN **
This is the primary detection of this object for this Frames ** Pipeline Run, and MAY be used in the primary survey, depending on the ** status of its other flags. **
DUPLICATE **
This object has a duplicate detection in an adjacent field ** (thus both objects lay in the overlap region between the two fields). ** A link will be made in the database to the duplicate detection. ** Such pairs of objects are useful for quality ** analysis internal to this run of the Frames Pipeline (i.e., the ** measured parameters should be all but exactly equal for each object ** in a pair). For all such pairs, exactly one object in the pair ** will also have the OK_RUN flag set. **
** Objects match if their separation is less than ** or equal to the specified match 'radius'. If more than one object in ** an adjacent field matches a given object, then only the closest match ** is merged. This returns a CHAIN of TSMATCHSTATs, with one set of ** statistics for each field in the run. The returned statistics give ** the row and column statistics in pixels, and the uncalibrated PSF ** magnitude statistics in magnitudes. Only GOOD objects classified STAR, ** GALAXY, or UNK are matched. **

** This verb should NOT be called directly. The verb "setObjectStatus" ** should be used; it calls this verb. ** ** ** **============================================================================ */ static char setObjectStatusDontCallCmd[] = "xSetObjectStatusDontCall"; static int setObjectStatusDontCallFlg = FTCL_ARGV_NO_LEFTOVERS; static ftclArgvInfo setObjectStatusDontCallTbl[] = { {NULL, FTCL_ARGV_HELP, NULL, NULL, "Set object status and self merge a Frames Pipeline Run\n"}, {"", FTCL_ARGV_STRING, NULL, NULL, "Handle of Frames Pipeline run"}, {"", FTCL_ARGV_DOUBLE, NULL, NULL, "Matching 'radius' (pixels)"}, {"", FTCL_ARGV_DOUBLE, NULL, NULL, "Compile matching statistics for nSigma detections (PSF, in r') only"}, {NULL, FTCL_ARGV_END, NULL, NULL, NULL} }; static int tclSetObjectStatusDontCall (ClientData clientData, Tcl_Interp *interp, int argc, char **argv) { ooHandle(arFramesPipeRun) *fp = NULL; char *framesName = NULL; double matchRadius, nSigma; int status; // Parse the arguments setObjectStatusDontCallTbl[1].dst = &framesName; setObjectStatusDontCallTbl[2].dst = &matchRadius; setObjectStatusDontCallTbl[3].dst = &nSigma; if ((status = shTclParseArgv(interp, &argc, argv, setObjectStatusDontCallTbl, setObjectStatusDontCallFlg, setObjectStatusDontCallCmd)) != FTCL_ARGV_SUCCESS) return status; // Get handle to frames pipeline run if (shTclAddrGetFromName (interp, framesName, (void **)&fp, "ooHandle(arFramesPipeRun)") != TCL_OK) { Tcl_SetResult(interp, "setObjectStatusDontCall: error getting handle to frames pipeline run", TCL_VOLATILE); return TCL_ERROR; } // Verify that this run has not previously been through setObjectStatus if ((*fp)->setObjectStatus().isValid() == oocTrue) { Tcl_SetResult(interp, "setObjectStatusDontCall: run was previously processed by setObjectStatus", TCL_VOLATILE); return TCL_ERROR; } // Get number of rows per frame. int camRow = 1; // Any photometric CCD row will do int rowsPerField = (*fp)->imagingRun()->program()->ccd(camRow,(*fp)->camCol())->nRows(); // HARDWIRED: number of overlap rows per frame. int overlapRows = 128; // Fetch the range of fields uint16 field0 = (*fp)->field0(); uint16 nFields = (*fp)->nFields(); // Allocate a CHAIN to return the matching statistics for each field. CHAIN *ch = shChainNewBlock("TSMATCHSTATS", nFields); if (ch == NULL) { Tcl_SetResult(interp, "setObjectStatusDontCall: couldn't allocate stats chain", TCL_VOLATILE); return TCL_ERROR; } // Allocate pointers for some dynamic structures TSMATCHSTATS *statsNew = NULL; TSMATCHSTATS *statsOld = NULL; ooTVArray(ooHandle(arObject)) *handlesNew = NULL; ooTVArray(ooHandle(arObject)) *handlesOld = NULL; TSOBJ *arrayNew = NULL; TSOBJ *arrayOld = NULL; // For each field ... CURSOR_T cursor = shChainCursorNew(ch); for (int32 field = field0; field < field0 + nFields; field++) { printf("%4d", field); fflush(stdout); // Save the pointers for the previous field if (field != field0) { arrayOld = arrayNew; handlesOld = handlesNew; statsOld = statsNew; } // Get the stats structure for this field if ((statsNew = (TSMATCHSTATS*) shChainWalk(ch, cursor, NEXT)) == NULL) { // Error. Clean up and return an error. shChainCursorDel(ch, cursor); taGenericChainDestroy(ch); if (handlesOld != NULL) delete handlesOld; if (arrayOld != NULL) shFree(arrayOld); char fieldString[10]; sprintf(fieldString, "%d", field); Tcl_AppendResult(interp, setObjectStatusDontCallCmd, ": error getting stats in field ", fieldString, NULL); return TCL_ERROR; } tsMatchStatsInit(statsNew); statsNew->field = field; // Read in the objects for this field, setting the SET, GOOD, and OK_RUN // flags as appropriate. Returns a pair of arrays of GOOD // objects in the field which MAY match objects in adjacent fields. // Each element in the "arrayNew" array contains an index into the // "handlesNew" array plus the object position in column and "scan row" // coordinates (i.e., row increases throughout the run). if (tsObjectStatusSet((*fp)->objectList(field), rowsPerField, overlapRows, matchRadius, nSigma, &handlesNew, &arrayNew, statsNew) == -1) { // Error. Clean up and return an error. shChainCursorDel(ch, cursor); taGenericChainDestroy(ch); if (handlesOld != NULL) delete handlesOld; if (arrayOld != NULL) shFree(arrayOld); shTclInterpAppendWithErrStack(interp); return TCL_ERROR; } // Skip remainder if this is the first field if (field == field0) continue; // Match objects between this and the previous field. Just search in // the lower overlap region +- the match radius. Sets the DUPLICATE // flag for all matches and guarantees that exactly one object in each // matched pair has the OK_RUN flag set. tsSelfMerge(handlesNew, arrayNew, handlesOld, arrayOld, matchRadius, field * rowsPerField - matchRadius, field * rowsPerField + overlapRows + matchRadius, statsNew); // Count the number of GOOD unmatched objects in the overlap regions for // the previous field (now that we've finished matching it). int size = handlesOld->size(); for (int k = 0; k < size; k++) { int status = (*handlesOld)[arrayOld[k].index]->status(0); if (! (status & AR_OBJECT_STATUS_OK_RUN) && ! (status & AR_OBJECT_STATUS_DUPLICATE)) { statsOld->nOverlap++; if (arrayOld[k].bright) statsOld->nOverlapBright++; } } // Delete info for previous field delete handlesOld; if (arrayOld != NULL) shFree(arrayOld); } shChainCursorDel(ch, cursor); printf("\n"); fflush(stdout); // Count the number of GOOD unmatched objects in the overlap regions for // the last field (now that we've finished matching it). int size = handlesNew->size(); for (int k = 0; k < size; k++) { int status = (*handlesNew)[arrayNew[k].index]->status(0); if (! (status & AR_OBJECT_STATUS_OK_RUN) && ! (status & AR_OBJECT_STATUS_DUPLICATE)) { statsNew->nOverlap++; if(arrayNew[k].bright) statsNew->nOverlapBright++; } } // Delete info for last field delete handlesNew; if (arrayNew != NULL) shFree(arrayNew); // Return a handle to the matching statistics chain char vName[HANDLE_NAMELEN]; if (shTclHandleNew(interp, vName, "CHAIN", (void *)ch) != TCL_OK) { // Error. Delete the chain and return an error. taGenericChainDestroy(ch); Tcl_SetResult(interp, "setObjectStatusDontCall: error allocating chain TCL handle", TCL_VOLATILE); return TCL_ERROR; } Tcl_SetResult(interp, vName, TCL_VOLATILE); return TCL_OK; } /*============================================================================ ** ** ** TCL VERB: resetObjectStatus ** ** ** Reset the SET, GOOD, and OK_RUN bits in the STATUS bit mask for ** each object in a frames pipeline run. This assumes that those bits were ** previously set, and that the new algorithm will only add new GOOD ** objects (that is, a previously GOOD object will never be declared ** BAD by the new algorithm). No border resolution between adjacent ** fields will be performed for the new OK_RUN objects. **

** This returns the number of new GOOD objects, new OK_RUN objects, old ** OK_RUN objects, and new OK_RUN objects in the border region between ** fields that could have possibly been matched to OK_RUN objects in ** adjacent fields but weren't and thus represents possible problems. ** ** ** **============================================================================ */ static char resetObjectStatusCmd[] = "resetObjectStatus"; static int resetObjectStatusFlg = FTCL_ARGV_NO_LEFTOVERS; static ftclArgvInfo resetObjectStatusTbl[] = { {NULL, FTCL_ARGV_HELP, NULL, NULL, "Reset object status for a Frames Pipeline Run\n"}, {"", FTCL_ARGV_STRING, NULL, NULL, "Handle of Frames Pipeline run"}, {"", FTCL_ARGV_INT, NULL, NULL, "Version of sky to reset (0=drilling, 1=current)"}, {NULL, FTCL_ARGV_END, NULL, NULL, NULL} }; static int tclresetObjectStatus (ClientData clientData, Tcl_Interp *interp, int argc, char **argv) { ooHandle(arFramesPipeRun) *fp = NULL; char *framesName = NULL; int skyVersion; int status; // Parse the arguments resetObjectStatusTbl[1].dst = &framesName; resetObjectStatusTbl[2].dst = &skyVersion; if ((status = shTclParseArgv(interp, &argc, argv, resetObjectStatusTbl, resetObjectStatusFlg, resetObjectStatusCmd)) != FTCL_ARGV_SUCCESS) return status; if (skyVersion != 0 && skyVersion != 1) { Tcl_SetResult(interp, "resetObjectStatus: 'skyVersion' must be 0 (drilling) or 1 (current)", TCL_VOLATILE); return TCL_ERROR; } // Get handle to frames pipeline run if (shTclAddrGetFromName (interp, framesName, (void **)&fp, "ooHandle(arFramesPipeRun)") != TCL_OK) { Tcl_SetResult(interp, "resetObjectStatus: error getting handle to frames pipeline run", TCL_VOLATILE); return TCL_ERROR; } // Verify that this run has previously been through setObjectStatus ooHandle(arProduct) setReport = (*fp)->setObjectStatus(); if (setReport.isValid() == oocFalse) { Tcl_SetResult(interp, "resetObjectStatus: run wasn't previously processed by setObjectStatus", TCL_VOLATILE); return TCL_ERROR; } // Get number of rows per frame. int camRow = 1; // Any photometric CCD row will do int rowsPerField = (*fp)->imagingRun()->program()->ccd(camRow,(*fp)->camCol())->nRows(); // HARDWIRED: number of overlap rows per frame. int overlapRows = 128; // Get the matchRadius used when this run was first processed // through setObjectStatus. double matchRadius = atof(setReport->parameter("matchRadius")); // As a sanity test, demand that matchRadius be greater than 0 arcsecs, // but less than (arbitrarily) 5 arcsecs. if (matchRadius <= 0. || matchRadius > 5.) { shErrStackPush("resetObjectStatus: bad value for merge parameter 'matchRadius'"); return SH_GENERIC_ERROR; } // Fetch the range of fields uint16 field0 = (*fp)->field0(); uint16 nFields = (*fp)->nFields(); // Initialize counts. int nNewGood = 0, nNewOKRun = 0, nOldOKRun = 0, nProblems = 0; // Issue warning printf("DANGER: YOU SHOULD NOT BE RUNNING THIS COMMAND!\n"); // For each field ... for (int32 field = field0; field < field0 + nFields; field++) { printf("%4d", field); fflush(stdout); // Reset objects for this field. if (tsObjectStatusReset((*fp)->objectList(field), rowsPerField, overlapRows, matchRadius, skyVersion, &nNewGood, &nNewOKRun, &nOldOKRun, &nProblems) == SH_GENERIC_ERROR) { // Error. Clean up and return an error. shTclInterpAppendWithErrStack(interp); return TCL_ERROR; } } printf("\n"); fflush(stdout); // Return a handle to the matching statistics chain char result[200]; sprintf(result, "{nNewGood %d} {nNewOKRun %d} {nOldOKRun %d} {nProblems %d}", nNewGood, nNewOKRun, nOldOKRun, nProblems); Tcl_AppendResult(interp, result, NULL); return TCL_OK; } /*============================================================================ ** ** ** TCL VERB: unsetObjectStatus ** ** ** Undo a "setObjectStatus" operation for a single Frames Pipeline Run. ** This clears the "arObjectStatus" flags for all objects in the ** run. It also removes all links for all objects in ** the run to other objects belonging to the run. It does not remove ** links to objects in other Frames Pipeline Runs. ** ** ** **============================================================================ */ static char runUnsetObjectStatusCmd[] = "unsetObjectStatus"; static int runUnsetObjectStatusFlg = FTCL_ARGV_NO_LEFTOVERS; static ftclArgvInfo runUnsetObjectStatusTbl[] = { {NULL, FTCL_ARGV_HELP, NULL, NULL, "Reset object status to UNDETERMINED for all objects in a Frames Pipeline Run, and break links to other objects in the same run\n"}, {"", FTCL_ARGV_STRING, NULL, NULL, "Handle of Frames Pipeline run"}, {"-force", FTCL_ARGV_CONSTANT, (void *) 1, NULL, "Force operation, even if run hasn't been through 'setObjectStatus'"}, {NULL, FTCL_ARGV_END, NULL, NULL, NULL} }; static int tclRunUnsetObjectStatus (ClientData clientData, Tcl_Interp *interp, int argc, char **argv) { ooHandle(arFramesPipeRun) *framesH = NULL; char *name = NULL; int force = 0; int status; // Parse the arguments runUnsetObjectStatusTbl[1].dst = &name; runUnsetObjectStatusTbl[2].dst = &force; if ((status = shTclParseArgv(interp, &argc, argv, runUnsetObjectStatusTbl, runUnsetObjectStatusFlg, runUnsetObjectStatusCmd)) != FTCL_ARGV_SUCCESS) return status; // Get handle to frames pipeline run if (shTclAddrGetFromName (interp, name, (void **)&framesH, "ooHandle(arFramesPipeRun)") != TCL_OK) { Tcl_AppendResult(interp, runUnsetObjectStatusCmd, ": error getting handle to run", NULL); return TCL_ERROR; } // Illegal to do so for a run which has already been merged against // other runs. if ((*framesH)->merge().isValid() == oocTrue) { Tcl_AppendResult(interp, runUnsetObjectStatusCmd, ": illegal: run has already been merged", NULL); return TCL_ERROR; } // Verify that this run has been through "setObjectStatus" if (! force) { if ((*framesH)->setObjectStatus().isValid() == oocFalse) { Tcl_AppendResult(interp, runUnsetObjectStatusCmd, ": run has not been through 'setObjectStatus'", NULL); return TCL_ERROR; } } // Fetch the range of fields uint16 field0 = (*framesH)->field0(); uint16 nFields = (*framesH)->nFields(); // For each field ... for (int32 field = field0; field < field0 + nFields; field++) { printf("%4d", field); fflush(stdout); // Set the status of all objects in the field to undetermined, and break // links to other objects in the same run. ooHandle(arObjectList) olist = (*framesH)->objectList(field); if (tsFieldUnsetObjectStatus(olist)!=SH_SUCCESS) { // Error. char fieldString[10]; sprintf(fieldString, "%d", field); Tcl_AppendResult(interp, runUnsetObjectStatusCmd, ": error in field ", fieldString, NULL); return TCL_ERROR; } } printf("\n"); fflush(stdout); // Remove link to "setObjectStatus" report ooHandle(arProduct) report = (*framesH)->setObjectStatus(); (*framesH)->delSetObjectStatus(); if (report.isValid() == oocTrue) { ooDelete(report); } // Return. return TCL_OK; } /*============================================================================ ** ** ** TCL VERB: selfMergeQA ** ** ** Generate matching statistics for when "setObjectStatus" was run ** on the specified Frames Pipeline Run. This returns a CHAIN of ** TSMATCHSTATS, with one set of statistics for each field in the run. ** The returned statistics give the row and column statistics in pixels, ** and the uncalibrated magnitude statistics (PSF for stars, petrosian ** for galaxies) in magnitudes. Only GOOD objects of class STAR, GALAXY, ** or UNK are counted. ** ** ** **============================================================================ */ static char selfMergeQACmd[] = "selfMergeQA"; static int selfMergeQAFlg = FTCL_ARGV_NO_LEFTOVERS; static ftclArgvInfo selfMergeQATbl[] = { {NULL, FTCL_ARGV_HELP, NULL, NULL, "Generate matching statistics between adjacent fields for a Frames Pipeline run\n"}, {"", FTCL_ARGV_STRING, NULL, NULL, "Handle of Frames Pipeline run"}, {"-firstField", FTCL_ARGV_INT, NULL, NULL, "First field"}, {"-lastField", FTCL_ARGV_INT, NULL, NULL, "Last field"}, {"-nSigma", FTCL_ARGV_DOUBLE, NULL, NULL, "Compile matching statistics for nSigma detections (PSF, in r') only (defaults to 10)"}, {"-noBlends", FTCL_ARGV_CONSTANT, (void *) 1, NULL, "Don't include blended objects (DFLAG_BLENDED || DFLAG_CHILD) in matching statistics and chains"}, {NULL, FTCL_ARGV_END, NULL, NULL, NULL} }; static int tclSelfMergeQA (ClientData clientData, Tcl_Interp *interp, int argc, char **argv) { ooHandle(arFramesPipeRun) *fp = NULL; char *name = NULL; int firstField = -1, lastField = -1; double nSigma = 10; int noBlends = 0; int status; // Parse the arguments selfMergeQATbl[1].dst = &name; selfMergeQATbl[2].dst = &firstField; selfMergeQATbl[3].dst = &lastField; selfMergeQATbl[4].dst = &nSigma; selfMergeQATbl[5].dst = &noBlends; if ((status = shTclParseArgv(interp, &argc, argv, selfMergeQATbl, selfMergeQAFlg, selfMergeQACmd)) != FTCL_ARGV_SUCCESS) return status; // Get handle to frames pipeline run if (shTclAddrGetFromName (interp, name, (void **)&fp, "ooHandle(arFramesPipeRun)") != TCL_OK) { Tcl_AppendResult(interp, selfMergeQACmd, ": error getting handle to run", NULL); return TCL_ERROR; } // Verify that this run has been through "setObjectStatus" if ((*fp)->setObjectStatus().isValid() == oocFalse) { Tcl_AppendResult(interp, selfMergeQACmd, ": run has not been through 'setObjectStatus'", NULL); return TCL_ERROR; } // Get number of rows per field. int camRow = 1; // Any photometric CCD row will do int rowsPerField = (*fp)->imagingRun()->program()->ccd(camRow,(*fp)->camCol())->nRows(); // Get/verify field range int field0 = (*fp)->field0(); int field1 = field0 + (*fp)->nFields() - 1; if (firstField > -1) { if (firstField < field0 || firstField > field1) { Tcl_AppendResult(interp, selfMergeQACmd, ": 'firstField' outside range of fields for this run", NULL); return TCL_ERROR; } } else { firstField = field0; } if (lastField > -1) { if (lastField < field0 || lastField > field1) { Tcl_AppendResult(interp, selfMergeQACmd, ": 'lastField' outside range of fields for this run", NULL); return TCL_ERROR; } } else { lastField = field1; } int nFields = lastField - firstField + 1; // Allocate a CHAIN to return the matching statistics for each field. CHAIN *ch = shChainNewBlock("TSMATCHSTATS", nFields); // For each field ... TSMATCHSTATS *stats; CURSOR_T cursor = shChainCursorNew(ch); for (int32 field = firstField; field <= lastField; field++) { printf("%4d", field); fflush(stdout); // Get the stats structure for this field if ((stats = (TSMATCHSTATS*) shChainWalk(ch, cursor, NEXT)) == NULL) { // Error. Clean up and return an error. shChainCursorDel(ch, cursor); taGenericChainDestroy(ch); char fieldString[10]; sprintf(fieldString, "%d", field); Tcl_AppendResult(interp, selfMergeQACmd, ": error getting stats in field ", fieldString, NULL); return TCL_ERROR; } stats->field = field; // Calculate the matching statistics for this field. Set the previous // field for the first field to itself --- this allows the subroutine to // work, but it will find no matches. if (tsFieldMatchStatsByPixel((*fp)->objectList(field), (*fp)->objectList(field == field0 ? field0 : field-1), rowsPerField, nSigma, stats, noBlends) != SH_SUCCESS) { // Error. Delete the chain and return an error. shChainCursorDel(ch, cursor); taGenericChainDestroy(ch); char fieldString[10]; sprintf(fieldString, "%d", field); Tcl_AppendResult(interp, selfMergeQACmd, ": error in field ", fieldString, NULL); return TCL_ERROR; } } shChainCursorDel(ch, cursor); printf("\n"); fflush(stdout); // Return a handle to the matching statistics chain char vName[HANDLE_NAMELEN]; if (shTclHandleNew(interp, vName, "CHAIN", (void *)ch) != TCL_OK) { // Error. Delete the chain and return an error. taGenericChainDestroy(ch); Tcl_AppendResult(interp, selfMergeQACmd, ": error allocating chain TCL handle", NULL); return TCL_ERROR; } Tcl_SetResult(interp, vName, TCL_VOLATILE); return TCL_OK; } /*============================================================================ ** ** ** TCL VERB: fieldPlot ** ** ** Plot the objects for a range of fields in a single Frames Pipeline Run. ** Two styles of plots are supported, the first to examine the results of ** "setObjectStatus", and the second the result of "resolve". ** In both plot styles, even numbered fields' objects are plotted in green, ** odd numbered fields in red. ** The symbol used for each point depends on the object's status. ** In the first style of plot (the default, useful for examine the results ** of "setObjectStatus"), ** triangle = matching object between fields (only one symbol is plotted ** for the pair, plotted in color of field to which the OK_RUN object in ** the pair belongs); dot = OK_RUN; plus = GOOD but not OK_RUN; ** cross = BAD. In the second style of plot, selected by specifying the ** "-resolve" flag (useful for examining the results of "resolve"), ** plus = unmatched primary, asterik = primary which matches one or more ** secondaries in another run, star = primary which matches one or more ** primaries in another run (which should be an error), ** triangle = secondary which matches a ** primary in another run, cross = secondary which does not match a ** primary. Primaries which have also been selected as spectroscopic ** targets are also marked with a circle. **

** In both types of plots, field boundaries are plotted in white. For the ** "-resolve" plots, scanline boundaries are plotted in yellow and stripe ** boundaries in blue. ** ** ** **============================================================================ */ static char fieldPlotCmd[] = "fieldPlot"; static int fieldPlotFlg = FTCL_ARGV_NO_LEFTOVERS; static ftclArgvInfo fieldPlotTbl[] = { {NULL, FTCL_ARGV_HELP, NULL, NULL, "Plot objects for a single Frames Pipeline Run\n"}, {"", FTCL_ARGV_STRING, NULL, NULL, "Handle of Frames Pipeline run"}, {"", FTCL_ARGV_INT, NULL, NULL, "0=drilling version, 1=current version"}, {"[device]", FTCL_ARGV_STRING, "/XWINDOW", NULL, "plotting device"}, {"-firstField", FTCL_ARGV_INT, NULL, NULL, "First field to plot (defaults to first field in the run)"}, {"-lastField", FTCL_ARGV_INT, NULL, NULL, "Last field to plot (defaults to last field in the run)"}, {"-resolve", FTCL_ARGV_CONSTANT, (void *) 1, NULL, "Plot symbols for debugging resolve"}, {NULL, FTCL_ARGV_END, NULL, NULL, NULL} }; static int tclFieldPlot (ClientData clientData, Tcl_Interp *interp, int argc, char **argv) { ooHandle(arFramesPipeRun) *fp = NULL; int skyVersion = 0; char *name = NULL; char *device = NULL; int firstField = -1, lastField = -1; int resolve = 0; int status; // Parse the arguments fieldPlotTbl[1].dst = &name; fieldPlotTbl[2].dst = &skyVersion; fieldPlotTbl[3].dst = &device; fieldPlotTbl[4].dst = &firstField; fieldPlotTbl[5].dst = &lastField; fieldPlotTbl[6].dst = &resolve; if ((status = shTclParseArgv(interp, &argc, argv, fieldPlotTbl, fieldPlotFlg, fieldPlotCmd)) != FTCL_ARGV_SUCCESS) return status; // Get handle to frames pipeline run if (shTclAddrGetFromName (interp, name, (void **)&fp, "ooHandle(arFramesPipeRun)") != TCL_OK) { Tcl_AppendResult(interp, fieldPlotCmd, ": error getting handle to run", NULL); return TCL_ERROR; } // Verify that this run has been through "setObjectStatus" if ((*fp)->setObjectStatus().isValid() == oocFalse) { Tcl_AppendResult(interp, fieldPlotCmd, ": run has not been through 'setObjectStatus'", NULL); return TCL_ERROR; } // If testing resolve, get astrometric calibration ooHandle(arAstromCalib) astrom; int ccd; double etaMin, etaMax, nuMin, nuMax, stripeNode, stripeIncl; if (resolve) { astrom = (*fp)->mergeAstrometricCalibration(); if (! astrom.isValid()) astrom = (*fp)->astrometricCalibration(); // All Frames Pipeline Runs are required when stuffed to have used r' as // their reference filter. Thus, we always want the TRANS structure for // the r' CCD. int camRow = 1; int camCol = (*fp)->camCol(); ccd = 10 * camRow + camCol; // Get the minimum and maximum survey latitude (eta) bounding this // stripe, and the minimum and maximum great circle latitude (nu) // bounding this scanline. int stripe = (*fp)->imagingRun()->stripe(); double etaCenter = at_etaMin + stripe * at_stripeSeparation; etaMin = etaCenter - at_stripeSeparation / 2.; etaMax = etaCenter + at_stripeSeparation / 2.; int scanline = (camCol - 1) * 2 + ((*fp)->imagingRun()->strip() == AR_STRIP_NORTH ? 1 : 0); nuMin = (scanline - 6) * at_scanSeparation; nuMax = nuMin + at_scanSeparation; stripeNode = 95.; stripeIncl = at_etaMin + stripe * at_stripeSeparation + at_surveyCenterDec; int southern = (stripe > 45); if (southern) stripeIncl -= 180.; } // Get/verify field range int field0 = (*fp)->field0(); int field1 = field0 + (*fp)->nFields() - 1; if (firstField > -1) { if (firstField < field0 || firstField > field1) { Tcl_AppendResult(interp, fieldPlotCmd, ": 'firstField' outside range of fields for this run", NULL); return TCL_ERROR; } } else { firstField = field0; } if (lastField > -1) { if (lastField < field0 || lastField > field1) { Tcl_AppendResult(interp, fieldPlotCmd, ": 'lastField' outside range of fields for this run", NULL); return TCL_ERROR; } } else { lastField = field1; } int nFields = lastField - firstField + 1; // Get number of rows and columns per field. int camRow = 1; // Any photometric CCD row will do ooHandle(arCCDProgram) program = (*fp)->imagingRun()->program()->ccd(camRow, (*fp)->camCol()); int rowsPerField = program->nRows(); // HARDWIRED: int colsPerField = 2048; // HARDWIRED: number of overlap rows per frame. int overlapRows = 128; // Open PGPLOT window if (cpgbeg(0, device, 1, 1) != 1) { Tcl_AppendResult(interp, fieldPlotCmd, ": couldn't open plotting device", NULL); return TCL_ERROR; } cpgask(0); // Plot the axes and label cpgenv(0, rowsPerField * nFields + overlapRows, -100, colsPerField+100, 1, 0); char label[100]; sprintf(label, "Run %d, camCol %d, rerun %d, fields %d-%d", (*fp)->run(), (*fp)->camCol(), (*fp)->rerun(), firstField, lastField); cpglab("x", "y", label); // Plot first field boundary float boundary = overlapRows/2.; cpgsci(1); cpgmove(boundary, 0); cpgdraw(boundary, colsPerField); // Plot points int offset = 0; cpgsfs(2); // plot symbols as outlines for (int field = firstField; field <= lastField; field++) { // Plot field boundary boundary = (field-firstField+1) * rowsPerField + overlapRows/2.; cpgsci(1); cpgmove(boundary, 0); cpgdraw(boundary, colsPerField); // If plotting for resolve, plot scanline and stripe limits double muMinOrig; TRANS trans; if (resolve) { // Plot scanline limits cpgsci(7); // yellow trans = astrom->astrotoolsTrans(ccd, field); double muMin = 0, muMax = 0, nu = 0; atTransApply(&trans, 'r', 0., 0., colsPerField / 2., 0., NULL, NULL, &muMin, NULL, &nu, NULL); atTransApply(&trans, 'r', rowsPerField + overlapRows, 0., colsPerField / 2., 0., NULL, NULL, &muMax, NULL, &nu, NULL); muMinOrig = muMin; // Parked scans don't scan the stripe great circle, so they // need to be converted. if (fabs(astrom->node() - stripeNode) > 0.00001 || fabs(astrom->incl() - stripeIncl) > 0.00001) { double ra, dec; atGCToEq(muMin, 0, &ra, &dec, astrom->node(), astrom->incl()); atEqToGC(ra, dec, &muMin, &nu, stripeNode, stripeIncl); atGCToEq(muMax, 0, &ra, &dec, astrom->node(), astrom->incl()); atEqToGC(ra, dec, &muMax, &nu, stripeNode, stripeIncl); } int nSteps = 100; double dMu = (muMax - muMin) / nSteps; muMin -= dMu; muMax += dMu; dMu = (muMax - muMin) / nSteps; for (int i = 0; i <= nSteps; i++) { double mu = muMin + i * dMu; nu = nuMin; // Parked scans don't scan the stripe great circle, so they // need to be converted. if (fabs(astrom->node() - stripeNode) > 0.00001 || fabs(astrom->incl() - stripeIncl) > 0.00001) { double ra, dec; atGCToEq(mu, nu, &ra, &dec, stripeNode, stripeIncl); atEqToGC(ra, dec, &mu, &nu, astrom->node(), astrom->incl()); atBound(&mu, muMinOrig, muMinOrig + 360.); } double row, col; atTransInverseApply(&trans, 'r', mu, 0., nu, 0., NULL, NULL, &row, NULL, &col, NULL); if (i == 0) cpgmove(row + offset, col); else cpgdraw(row + offset, col); } for (i = 0; i <= nSteps; i++) { double mu = muMin + i * dMu; nu = nuMax; // Parked scans don't scan the stripe great circle, so they // need to be converted. if (fabs(astrom->node() - stripeNode) > 0.00001 || fabs(astrom->incl() - stripeIncl) > 0.00001) { double ra, dec; atGCToEq(mu, nu, &ra, &dec, stripeNode, stripeIncl); atEqToGC(ra, dec, &mu, &nu, astrom->node(), astrom->incl()); atBound(&mu, muMinOrig, muMinOrig+360.); } double row, col; atTransInverseApply(&trans, 'r', mu, 0., nu, 0., NULL, NULL, &row, NULL, &col, NULL); if (i == 0) cpgmove(row + offset, col); else cpgdraw(row + offset, col); } // Plot stripe limits cpgsci(4); // blue double lambdaMin, lambdaMax, eta; atGCToSurvey(muMin, 0, stripeNode, stripeIncl, &lambdaMin, &eta); atGCToSurvey(muMax, 0, stripeNode, stripeIncl, &lambdaMax, &eta); double dLambda = (lambdaMax - lambdaMin) / nSteps; for (i = 0; i <= nSteps; i++) { double lambda = lambdaMin + i * dLambda; double mu, nu; atSurveyToGC(lambda, etaMin, astrom->node(), astrom->incl(), &mu, &nu); atBound(&mu, muMinOrig, muMinOrig+360.); double row, col; atTransInverseApply(&trans, 'r', mu, 0., nu, 0., NULL, NULL, &row, NULL, &col, NULL); if (i == 0) cpgmove(row + offset, col); else cpgdraw(row + offset, col); } for (i = 0; i <= nSteps; i++) { double lambda = lambdaMin + i * dLambda; double mu, nu; atSurveyToGC(lambda, etaMax, astrom->node(), astrom->incl(), &mu, &nu); atBound(&mu, muMinOrig, muMinOrig+360.); double row, col; atTransInverseApply(&trans, 'r', mu, 0., nu, 0., NULL, NULL, &row, NULL, &col, NULL); if (i == 0) cpgmove(row + offset, col); else cpgdraw(row + offset, col); } } // Get object list ooHandle(arObjectList) list = (*fp)->objectList(field); // Plot segment boundary for each segments for which this list is the // first or last field in the segment. if (resolve) { ooItr(arSegment) segments; char predicate[100]; sprintf(predicate, "chunk_.skyVersion_ == %d", skyVersion); list->segments(segments, predicate); while (segments.next()) { if (field == segments->field0() || field == segments->field0() + segments->nFields() - 1) { cpgsci(7); double mu; if (field == segments->field0()) mu = segments->startMu() * at_asec2Deg; else mu = segments->endMu() * at_asec2Deg; int nSteps = 100; double dNu = (nuMax - nuMin) / nSteps; for (int i = 0; i <= nSteps; i++) { double mu1 = mu; double nu1 = nuMin + i * dNu; if (fabs(astrom->node() - stripeNode) > 0.00001 || fabs(astrom->incl() - stripeIncl) > 0.00001) { double ra, dec; atGCToEq(mu1, nu1, &ra, &dec, stripeNode,stripeIncl); atEqToGC(ra, dec, &mu1, &nu1, astrom->node(), astrom->incl()); atBound(&mu1, muMinOrig, muMinOrig+360.); } double row, col; atTransInverseApply(&trans, 'r', mu1, 0., nu1, 0., NULL, NULL, &row, NULL, &col, NULL); if (i == 0) cpgmove(row + offset, col); else cpgdraw(row + offset, col); } } } } // Set color of objects in this and next field. int thisColor, nextColor; if ((field / 2) * 2 == field) { thisColor = 2; nextColor = 3; } else { thisColor = 3; nextColor = 2; } // Get iterator to objects for this field ooItr(arObject) itr; if (list->objects(itr, oocReadOnly) != oocSuccess) { Tcl_SetResult(interp, "selfMergePlot: couldn't initialize iterator over objects", TCL_VOLATILE); return TCL_ERROR; } // For each object ... while (itr.next()) { // Set color to default color for this field. cpgsci(thisColor); // Get row and col float x = itr->objc_rowc().value() + offset; float y = itr->objc_colc().value(); // Get object status int status = itr->status(0); int symbol; int secondSymbol = 0; if (resolve) { // Only plot OK_RUN objects if (! (status & AR_OBJECT_STATUS_OK_RUN)) continue; int resolveStatus = itr->status(skyVersion); // Plot unresolved objects as dots if (! (resolveStatus & AR_OBJECT_STATUS_RESOLVED)) symbol = 1; else { // Check for matching objects. int matchStatus = -1; ooItr(arObject) matches; itr->matches(matches); while (matches.next()) { // Must test because run 125, rerun 4, camCol 5 contains // some matches to non-existent objects. For example: // 125-4-400-710 -> 125-5-5-400-716 -> # 1929-392-4-44 // -> ??? -> # 2577-392-5-31 // -> 125-7-5-400-698 -> # 2623-392-5-23 // The database #2577 doesn't exist. if (! matches.isValid()) { printf("WARNING: invalid handle to a matched object found for object %d-%d-%d-%d-%d\n", itr->run(), itr->rerun(), itr->camCol(), itr->field(), itr->id()); continue; } if (matches->status(skyVersion) & AR_OBJECT_STATUS_PRIMARY) { matchStatus = 1; break; } else if (matches->status(skyVersion) & AR_OBJECT_STATUS_SECONDARY) matchStatus = 0; } // Plot unmatched primaries as a plus, primaries which // match a secondary as an asterik, primaries which match // a primary as a star, secondaries which match // a primary as a triangle, unmatched secondaries as a // cross, and objects which are neither primary nor // secondary as a dot. if (resolveStatus & AR_OBJECT_STATUS_PRIMARY) { if (matchStatus == 1) symbol = 12; else if (matchStatus == 0) symbol = 3; else symbol = 2; // If a spectroscopic target, mark with a circle if a // QSO, else a square if a galaxy, else a diamond. if (resolveStatus & AR_OBJECT_STATUS_TARGET && \ itr->priority(skyVersion) > 0) { uint32 primaryType = itr->primaryType(skyVersion); if (primaryType & (AR_TARGET_QSO_HIZ | AR_TARGET_QSO_CAP | AR_TARGET_QSO_SKIRT | AR_TARGET_QSO_FIRST_CAP | AR_TARGET_QSO_FIRST_SKIRT)) secondSymbol = 4; else if (primaryType & (AR_TARGET_GALAXY_RED | AR_TARGET_GALAXY | AR_TARGET_GALAXY_BIG | AR_TARGET_GALAXY_BRIGHT_CORE)) secondSymbol = 6; else secondSymbol = 11; } } else if (resolveStatus & AR_OBJECT_STATUS_SECONDARY) { if (matchStatus == 1) symbol = 7; else symbol = 5; } else { symbol = 1; } } } else { if (! (status & AR_OBJECT_STATUS_GOOD)) symbol = 5; // plot bad objects as crosses else if (status & AR_OBJECT_STATUS_DUPLICATE) { // Plot matching objects just once. Only plot when in // upper overlap region, except in first field. Plot as // triangles. if (itr->objc_rowc().value() < rowsPerField / 2. && field != firstField) continue; symbol = 7; if (! (status & AR_OBJECT_STATUS_OK_RUN)) cpgsci(nextColor); } else if (status & AR_OBJECT_STATUS_OK_RUN) symbol = 1; // plot OK_RUN objects as dots else symbol = 2; // plot non-matching !OK_RUN objects as pluses } // Plot the point cpgpt(1, &x, &y, symbol); if (secondSymbol > 0) cpgpt(1, &x, &y, secondSymbol); } offset = offset + rowsPerField; } cpgask(1); // Return nothing Tcl_SetResult(interp, "", TCL_VOLATILE); return TCL_OK; } /*============================================================================ ** ** ** TCL VERB: xRunMergeDontCall ** ** ** Merge the OK_RUN objects (only those classified as STAR, GALAXY, or ** UNK) in a Frames Pipeline Run with OK_RUN stars and galaxies in all ** other previously merged Frames Pipeline Runs. Objects match if their ** separation is less than or ** equal to the specified match 'radius'. If an object in a single ** previously-merged Frames Pipeline Run matches more than one object in ** the target run, then only the closest match is merged (an object in ** the target run can match more than one object in each of the ** previously-merged run). **

** The specified astrometric calibration is used to derive calibrated ** positions for matching. ** DCR corrections are applied using cosmic colors, since ** a final photometric calibration can not be guaranteed to exist for ** the run at this time. **

** This verb should NOT be called directly. The verb "runMerge" should ** should be used; it calls this verb. ** ** ** **============================================================================ */ static char runMergeDontCallCmd[] = "xRunMergeDontCall"; static int runMergeDontCallFlg = FTCL_ARGV_NO_LEFTOVERS; static ftclArgvInfo runMergeDontCallTbl[] = { {NULL, FTCL_ARGV_HELP, NULL, NULL, "Merge the objects in a Frames Pipeline Run with objects in all other previously merged Frames Pipeline Runs\n"}, {"", FTCL_ARGV_STRING, NULL, NULL, "Handle of Frames Pipeline run (ooHandle(arFramesPipeRun))"}, {"", FTCL_ARGV_STRING, NULL, NULL, "Handle of astrometric calibration to use (ooHandle(arAstromCalib))"}, {"", FTCL_ARGV_DOUBLE, NULL, NULL, "Match 'radius' (arcsecs)"}, {NULL, FTCL_ARGV_END, NULL, NULL, NULL} }; static int tclRunMergeDontCall (ClientData clientData, Tcl_Interp *interp, int argc, char **argv) { ooHandle(arFramesPipeRun) *framesH = NULL; ooHandle(arAstromCalib) *astromH = NULL; char *framesName = NULL; char *astromName = NULL; double matchRadius; int status; // Parse the arguments runMergeDontCallTbl[1].dst = &framesName; runMergeDontCallTbl[2].dst = &astromName; runMergeDontCallTbl[3].dst = &matchRadius; if ((status = shTclParseArgv(interp, &argc, argv, runMergeDontCallTbl, runMergeDontCallFlg, runMergeDontCallCmd)) != FTCL_ARGV_SUCCESS) return status; // Get handle to frames pipeline run if (shTclAddrGetFromName (interp, framesName, (void **)&framesH, "ooHandle(arFramesPipeRun)") != TCL_OK) { Tcl_AppendResult(interp, runMergeDontCallCmd, ": error getting handle to frames pipeline run", NULL); return TCL_ERROR; } // Verify that this run has not previously been merged if ((*framesH)->merge().isValid() == oocTrue) { Tcl_AppendResult(interp, runMergeDontCallCmd, ": run '", framesName, "' has previously been merged", NULL); return TCL_ERROR; } // Verify that this run has been through setObjectStatus. if ((*framesH)->setObjectStatus().isValid() == oocFalse) { Tcl_AppendResult(interp, runMergeDontCallCmd, ": run '", framesName, "' has not been through setObjectStatus", NULL); return TCL_ERROR; } // Get handle to astrometric calibration if (shTclAddrGetFromName (interp, astromName, (void **)&astromH, "ooHandle(arAstromCalib)") != TCL_OK) { Tcl_AppendResult(interp, runMergeDontCallCmd, ": error getting handle to astrometric calibration", NULL); return TCL_ERROR; } // Verify that astrometric calibration is for the target imaging run and // spans the full range of fields if ((*astromH)->run() != (*framesH)->run()) { Tcl_AppendResult(interp, runMergeDontCallCmd, ": astrometric calibration for wrong imaging run", NULL); return TCL_ERROR; } int astromField0 = (*astromH)->field0(); int framesField0 = (*framesH)->field0(); if (astromField0 > framesField0 || astromField0+(*astromH)->nFields() < framesField0+(*framesH)->nFields()) { Tcl_AppendResult(interp, runMergeDontCallCmd, ": astrometric calibration doesn't cover field range of the frames pipeline run", NULL); return TCL_ERROR; } // Fetch the global TCL variable "photo_frames" for container with all // Frames Pipeline Runs. Then fetch the associated handle. char *photo_frames = Tcl_GetVar(interp, "photo_frames", TCL_GLOBAL_ONLY); if (photo_frames == NULL) { Tcl_AppendResult(interp, runMergeDontCallCmd, ": global variable $photo_frames doesn't exist", NULL); return TCL_ERROR; } ooHandle(ooContObj) *framesContH; if (shTclAddrGetFromName (interp, photo_frames, (void **)&framesContH, "ooHandle(ooContObj)") != TCL_OK) { Tcl_AppendResult(interp, runMergeDontCallCmd, ": error getting handle to $photo_frames", NULL); return TCL_ERROR; } // Fetch great circle parameters for target run double equinox = (*astromH)->equinox(); double node = (*astromH)->node(); double incl = (*astromH)->incl(); // Initialize iterator over all frames pipeline runs ooItr(arFramesPipeRun) itr; if (itr.scan(*framesContH, oocRead) != oocSuccess) { Tcl_AppendResult(interp, runMergeDontCallCmd, ": couldn't initalize iterator over frames pipeline runs", NULL); return TCL_ERROR; } // Set up cosmic color for use by atTransApply(). taTransCosmicSet(cosmicColor, cosmicScatter); // Declare some needed variables TSCALIBOBJ *array; ooTVArray(ooHandle(arObject)) *handles; int nObjects; double minNu, maxNu; // For each frames pipeline run ... int first = 1; while (itr.next()) { // Skip if this is the same run as the target run if (itr == *framesH) continue; // Skip if this run has not yet been merged itself. ooHandle(arProduct) report; report = itr->merge(report); if (report.isValid() == oocFalse) continue; // TODO: we should test whether this run could possibly overlap // the target run, and skip it if it can't. // If we haven't, fetch the array of sorted (by increasing mu) // calibrated objects for the target run. if (first == 1) { if ((nObjects = tsRunToCalibArray(*framesH, *astromH, &handles, &array, &minNu, &maxNu)) == -1) { shTclInterpAppendWithErrStack(interp); return TCL_ERROR; } if (nObjects == 0) { delete handles; shFree(array); Tcl_SetResult(interp, "", TCL_VOLATILE); return TCL_OK; } first = 0; } // Merge the target run against this run if (tsRunMergeOne(handles, array, nObjects, minNu, maxNu, node, incl, equinox, itr, matchRadius) != SH_SUCCESS) { delete handles; shFree(array); shTclInterpAppendWithErrStack(interp); return TCL_ERROR; } } // Delete the array if (first == 0) { delete handles; shFree(array); } // Set the calibrations used (*framesH)->setMergeAstrometricCalibration(*astromH); // Return Tcl_SetResult(interp, "", TCL_VOLATILE); return TCL_OK; } /*============================================================================ ** ** ** TCL VERB: runUnmerge ** ** ** Undo a "merge" operation of a single Frames Pipeline run. This ** unmerges a Frames Pipeline run from other merged Frames Pipeline runs. ** It breaks all links between objects in the Frames Pipeline run with ** objects in other Frames Pipeline runs. It does not remove ** links to objects in the same Frames Pipeline Run. ** ** ** **============================================================================ */ static char runUnmergeCmd[] = "runUnmerge"; static int runUnmergeFlg = FTCL_ARGV_NO_LEFTOVERS; static ftclArgvInfo runUnmergeTbl[] = { {NULL, FTCL_ARGV_HELP, NULL, NULL, "Remove links to objects in other Frames Pipeline runs\n"}, {"", FTCL_ARGV_STRING, NULL, NULL, "Handle of Frames Pipeline run"}, {"-force", FTCL_ARGV_CONSTANT, (void *) 1, NULL, "Force operation, even if run hasn't been through 'merge'"}, {NULL, FTCL_ARGV_END, NULL, NULL, NULL} }; static int tclRunUnmerge (ClientData clientData, Tcl_Interp *interp, int argc, char **argv) { ooHandle(arFramesPipeRun) *framesH = NULL; char *name = NULL; int force = 0; int status; // Parse the arguments runUnmergeTbl[1].dst = &name; runUnmergeTbl[2].dst = &force; if ((status = shTclParseArgv(interp, &argc, argv, runUnmergeTbl, runUnmergeFlg, runUnmergeCmd)) != FTCL_ARGV_SUCCESS) return status; // Get handle to frames pipeline run if (shTclAddrGetFromName (interp, name, (void **)&framesH, "ooHandle(arFramesPipeRun)") != TCL_OK) { Tcl_AppendResult(interp, runUnmergeCmd, ": error getting handle to run", NULL); return TCL_ERROR; } // Verify that this run does not belong to any resolved segments. ooItr(arSegment) itr; if ((*framesH)->segments(itr) == oocError) { Tcl_AppendResult(interp, runUnmergeCmd, ": error initialize iterator over segments", NULL); return TCL_ERROR; } if (itr.next() == oocTrue) { Tcl_AppendResult(interp, runUnmergeCmd, ": run belongs to one or more resolved segments", NULL); return TCL_ERROR; } // Verify that this run has been through "merge" ooHandle(arProduct) mergeReport = (*framesH)->merge(); if (! force) { if (mergeReport.isValid() == oocFalse) { Tcl_AppendResult(interp, runUnmergeCmd, ": run has not been merged", NULL); return TCL_ERROR; } } // Fetch the range of fields uint16 field0 = (*framesH)->field0(); uint16 nFields = (*framesH)->nFields(); // For each field ... for (int32 field = field0; field < field0 + nFields; field++) { printf("%4d", field); fflush(stdout); // Remove links to objects in other runs. ooHandle(arObjectList) olist = (*framesH)->objectList(field); if (tsFieldUnmerge(olist) != SH_SUCCESS) { // Error. char fieldString[10]; sprintf(fieldString, "%d", field); Tcl_AppendResult(interp, runUnmergeCmd, ": error in field ", fieldString, NULL); return TCL_ERROR; } } printf("\n"); fflush(stdout); // Remove link to "merge" report and the calibrations used for the merge, // and delete the merge report from the database. (*framesH)->delMerge(); (*framesH)->delMergeAstrometricCalibration(); if (mergeReport.isValid() == oocTrue) { ooDelete(mergeReport); } // Return. return TCL_OK; } /*============================================================================ ** ** ** TCL VERB: runMergeQA ** ** ** Generate matching statistics for the matching between two Frames ** Pipeline runs. This returns a CHAIN of TSMATCHSTAT, with one set of ** statistics for each field in the first run. If the astrometric or ** photometric calibrations to use for either run are not specified, then ** the "current" calibrations are used. ** Position statistics are in calibrated great circle coordinates ** (arcsecs, where the coordinate system of "run1" is used). **

** Only OK_RUN objects of class STAR, GALAXY, or UNK are used, with ** separate statistics calculated for stars and galaxies. Objects to ** consider are restricted to the specified magnitude range. PSF ** magnitudes are used for stars and unknowns, petrosian magnitudes for ** galaxies (this can be overridden by specifying a magnitude to be used ** regardless of image classification). Objects are further ** restricted to those with specified flags turned on and specified flags ** turned off (separate sets of flags for the two different runs); these ** flags can be either the objc_flags or the flags in the specified ** filter. Flags are specified as a TCL list of ascii values. For ** example, to examine only run 1 objects which are EDGE objects and not ** blended (neither BLEND nor CHILD flags set), matched up against EDGE ** objects in run 2: ** ** 

**          runMergeQA ... -flags1On EDGE -flags1Off {BLEND CHILD}
**                         -flags2On EDGE
**      

**	Be aware then if the flags are fetched from an individual filter,
**	the entire object must be read, rather than the lite object, thus
**	slowly execution time down considerably.
**	If the "-obj" flag is specified, then three chains containing each
**	matched pair are returned, one for star-star pairs, one for
**      galaxy-galaxy pairs, and one for all other pairs (mismatches and pairs
**      containing UNK objects).  In this
**	case, the procedure returns a TCL list, where the first element is
**	the stats chain, the second element is the matched stars chain,
**	the third element is the matched galaxies chain, and the fourth
**      element is the remaining pairs chain.  If the "-full" option is
**	specified, then two additional elements are returned; the fifth
**	element is the field info chain for run 1, and the sixth element is
**	the field info chain for run 2.
**	

**      If "missingDistance" is specified,
**	then bright OK_RUN objects of class STAR, GALAXY, or UNK in the first
**	run which don't match an object in the second run but should have are
**	returned on that chain.  "Should have" means that they should be
**	located at least "missingDistance" pixels from all boundaries for a
**	field in the other run.  The chain of missing objects will be appended
**	to the return list.
**	

**	If "[run2]" is not specified, then
**	no matching occurs.  Only the non-matching statistics are calculated,
**	with the sample defined as described above.  The statistics chains is
**	returned.  If "-full" is specified, then four additional chains are
**	returned (again as a TCL list, with the statistics chain as the first
**	element in a give element list).  The second, third, and fourth
**	chain elements are TA_CALIB_OBJ chains containing all objects meeting
**	the selection criterion of class STAR, GALAXY, and UNK, respectively,
**	and the fifth element is the chain of field info for the run.
**	

**	If "-raw" is specified, then lengths, magnitudes, surface brightnesses
**	and angles are not calibrated, in the full versions of objects only.
**	Celestial coordinates and extinction
**	values are still calculated, however.
**
**
**
**============================================================================
*/
static char runMergeQACmd[] = "runMergeQA";
static int runMergeQAFlg = FTCL_ARGV_NO_LEFTOVERS;
static ftclArgvInfo runMergeQATbl[] = {
   {NULL, FTCL_ARGV_HELP, NULL, NULL,
    "Generate matching statistics between two Frames Pipeline runs\n"},
   {"", FTCL_ARGV_STRING, NULL, NULL,
    "Handle of first Frames Pipeline run (ooHandle(arFramesPipeRun))"},
   {"[run2]", FTCL_ARGV_STRING, NULL, NULL,
    "Handle of second Frames Pipeline run (ooHandle(arFramesPipeRun))"},
   {"-astrom1", FTCL_ARGV_STRING, NULL, NULL,
    "Handle of astrometric calibration to use for run 1 (ooHandle(arAstromCalib))"},
   {"-astrom2", FTCL_ARGV_STRING, NULL, NULL,
    "Handle of astrometric calibration to use for run 2 (ooHandle(arAstromCalib))"},
   {"-photo1", FTCL_ARGV_STRING, NULL, NULL,
    "Handle of photometric calibration to use for run 1 (ooHandle(arFinalPhotoCalib))"},
   {"-photo2", FTCL_ARGV_STRING, NULL, NULL,
    "Handle of photometric calibration to use for run 2 (ooHandle(arFinalPhotoCalib))"},
   {"-firstField", FTCL_ARGV_INT, NULL, NULL, "First field"},
   {"-lastField", FTCL_ARGV_INT, NULL, NULL, "Last field"},
   {"-min", FTCL_ARGV_DOUBLE, "15", NULL,
    "Minimum magnitude for 'bright' objects"},
   {"-max", FTCL_ARGV_DOUBLE, "19", NULL,
    "Maximum magnitude for 'bright' objects"},
   {"-filter", FTCL_ARGV_STRING, "r", NULL,
    "Filter used for 'bright' magnitude test (ugriz)"},
   {"-psf", FTCL_ARGV_CONSTANT, (void *) 1, NULL,
    "Select bright objects based on PSF magnitudes, regardless of classification"},
   {"-petro", FTCL_ARGV_CONSTANT, (void *) 1, NULL,
    "Select bright objects based on petrosian magnitudes, regardless of classification"},
   {"-model", FTCL_ARGV_CONSTANT, (void *) 1, NULL,
    "Select bright objects based on model magnitudes, regardless of classification"},
   {"-obj", FTCL_ARGV_CONSTANT, (void *) 1, NULL,
    "Return chains for matched stars and galaxies"},
   {"-full", FTCL_ARGV_CONSTANT, (void *) 1, NULL,
    "If returning chains for matched stars and galaxies, return the full objects, rather than just mags and positions (forces -obj)"},
   {"-flagsFilter", FTCL_ARGV_STRING, NULL, NULL, "Filter to fetch flags from (u, g, r, i, or z --- defaults to objc_flags"},
   {"-flags1On", FTCL_ARGV_STRING, NULL, NULL,
    "Flags which must be on for objects in run 1"},
   {"-flags1Off", FTCL_ARGV_STRING, NULL, NULL,
    "Flags which must be off for objects in run 1"},
   {"-flags2On", FTCL_ARGV_STRING, NULL, NULL,
    "Flags which must be on for objects in run 2"},
   {"-flags2Off", FTCL_ARGV_STRING, NULL, NULL,
    "Flags which must be off for objects in run 2"},
   {"-raw", FTCL_ARGV_CONSTANT, (void *) 1, NULL,
    "Don't calibrate lengths, magnitudes, surface brightesses, and angles in full version of objects"},
   {"-missingDistance", FTCL_ARGV_DOUBLE, NULL, NULL,
    "Return chain of missing objects using this distance from the other runs fields boundaries to define missing (pixels)"},
   {NULL, FTCL_ARGV_END, NULL, NULL, NULL}
};

static int
tclRunMergeQA (ClientData clientData,
			 Tcl_Interp *interp,
			 int argc,
			 char **argv)
{
   /*adjacent filter to use in colorterms*/
   const int adj[nFilters]  = {1, 2, 3, 4, 3};
   /* sign to use for color terms */
   const int sign[nFilters] = {1, 1, 1, 1, -1}; 

   ooHandle(arFramesPipeRun) *fp[2] = {NULL, NULL};
   ooHandle(arAstromCalib) *astrom[2] = {NULL, NULL};
   ooHandle(arAstromCalib) astromH[2];
   ooHandle(arFinalPhotoCalib) *final[2] = {NULL, NULL};
   ooHandle(arFinalPhotoCalib) finalH[2];
   char *fpName[2] = {NULL, NULL};
   char *astromName[2] = {NULL, NULL};
   char *finalName[2] = {NULL, NULL};
   int firstField = -1, lastField = -1;
   double minMag, maxMag;
   char *filter;
   int psfMags = 0;
   int petroMags = 0;
   int modelMags = 0;
   int mags;
   double expTime[2];
   int objChains = 0;
   int full = 0;
   char *flagsFilter = " ";
   char *flagsString[2][2] = {NULL, NULL, NULL, NULL};
   int raw = 0;
   double missingDistance = -1;
   int status;

   // Parse the arguments
   runMergeQATbl[1].dst = &fpName[0];
   runMergeQATbl[2].dst = &fpName[1];
   runMergeQATbl[3].dst = &astromName[0];
   runMergeQATbl[4].dst = &astromName[1];
   runMergeQATbl[5].dst = &finalName[0];
   runMergeQATbl[6].dst = &finalName[1];
   runMergeQATbl[7].dst = &firstField;
   runMergeQATbl[8].dst = &lastField;
   runMergeQATbl[9].dst = &minMag;
   runMergeQATbl[10].dst = &maxMag;
   runMergeQATbl[11].dst = &filter;
   runMergeQATbl[12].dst = &psfMags;
   runMergeQATbl[13].dst = &petroMags;
   runMergeQATbl[14].dst = &modelMags;
   runMergeQATbl[15].dst = &objChains;
   runMergeQATbl[16].dst = &full;
   runMergeQATbl[17].dst = &flagsFilter;
   runMergeQATbl[18].dst = &(flagsString[0][0]);
   runMergeQATbl[19].dst = &(flagsString[0][1]);
   runMergeQATbl[20].dst = &(flagsString[1][0]);
   runMergeQATbl[21].dst = &(flagsString[1][1]);
   runMergeQATbl[22].dst = &raw;
   runMergeQATbl[23].dst = &missingDistance;
   if ((status = shTclParseArgv(interp, &argc, argv, runMergeQATbl,
				runMergeQAFlg, runMergeQACmd))
       != FTCL_ARGV_SUCCESS)
     return status;
   if (full) objChains = 1; // -full implies -obj
   if (psfMags + petroMags + modelMags > 1)
     {
       Tcl_AppendResult(interp, runMergeQACmd,
			": only one of -psfMags, -petroMags, and -modelMags can bespecified", NULL);
       return TCL_ERROR;
     }
   if (psfMags)
     mags = 1;
   else if (petroMags)
     mags = 2;
   else if (modelMags)
     mags = 3;
   else
     mags = 0;
   
   // If run2 not specified, then we're not matching.
   int match = 1;
   if (fpName[1] == NULL) match = 0;
   if (missingDistance >= 0 && match == 0)
     {
       Tcl_AppendResult(interp, runMergeQACmd,
			": can't specify 'missingDistance' when not matching",
			NULL);
       return TCL_ERROR;
     }

   // Decode flags strings
   GenericEnum* enum1 = GenericEnum::Find("AR_DETECTION_FLAG");
   shAssert(enum1 != NULL);
   GenericEnum* enum2 = GenericEnum::Find("AR_DETECTION_FLAG2");
   shAssert(enum2 != NULL);
   char *prefix1 = "AR_DFLAG_";
   char *prefix2 = "AR_DFLAG2_";
   int flags[2][2][2] = {0, 0, 0, 0, 0, 0}; // [run1/2][enum1/2][on/off]
   for (int i = 0; i < 2; i++)
     for (int j = 0; j < 2; j++)
       {
	 if (flagsString[i][j] == NULL) continue;
	 char **listArgv;
	 int nFlags = 0;
	 if (Tcl_SplitList(interp, flagsString[i][j], &nFlags, &listArgv)
	     != TCL_OK)
	   {
	     Tcl_AppendResult(interp, runMergeQACmd,
			      ": couldn't parse flags string", NULL);
	     return TCL_ERROR;
	   }
	 char buffer[100];
	 for (int k = 0; k < nFlags; k++)
	   {
	     if (strlen(listArgv[k]) > 40)
	       {
		 Tcl_AppendResult(interp, runMergeQACmd,
				  ": illegal flag '", listArgv[k], "'", NULL);
		 return TCL_ERROR;
	       }
	     sprintf(buffer, "%s%s", prefix1, listArgv[k]);
	     if (enum1->validKey(buffer) == TRUE)
	       flags[i][0][j] |= enum1->asInteger(buffer);
	     else
	       {
		 sprintf(buffer, "%s%s", prefix2, listArgv[k]);
		 if (enum2->validKey(buffer) == TRUE)
		   flags[i][1][j] |= enum2->asInteger(buffer);
		 else
		   {
		     Tcl_AppendResult(interp, runMergeQACmd,
				      ": illegal flag '", listArgv[k], "'",
				      NULL);
		     return TCL_ERROR;
		   }
	       }
	   }
       }

   // For each input frames pipeline run ...
   for (i = 0; i < 2; i++)
     {
       if (i == 1 && ! match) continue;

       // Get handle to the run
       if (shTclAddrGetFromName (interp, fpName[i], (void **)&fp[i],
				 "ooHandle(arFramesPipeRun)") != TCL_OK)
	 {
	   Tcl_AppendResult(interp, runMergeQACmd,
			    ": error getting handle to run '", fpName[i],
			    ",", NULL);
	   return TCL_ERROR;
	 }
   
       // Verify that this run has been merged
       if ((*fp[i])->merge().isValid() == oocFalse && match)
	 {
	   Tcl_AppendResult(interp, runMergeQACmd, ": run '", fpName[i],
			    "' has not been merged", NULL);
	   return TCL_ERROR;
	 }
       
       // Get handle to the astrometric calibration
       char buffer[10];
       sprintf(buffer, "-astrom%d", i+1);
       if (ftcl_ArgIsPresent(argc, argv, buffer, runMergeQATbl) == 1)
	 {
	   // If one was specified, use it.
	   if (shTclAddrGetFromName (interp, astromName[i],
				     (void **)&astrom[i],
				     "ooHandle(arAstromCalib)") != TCL_OK)
	     {
	       Tcl_AppendResult(interp, runMergeQACmd,
				": error getting handle to astrom run '",
				astromName[i], "'", NULL);
	       return TCL_ERROR;
	     }
	   if ((*astrom[i])->run() != (*fp[i])->run())
	     {
	       Tcl_AppendResult(interp, runMergeQACmd,
				": astrometric calibration '",astromName[i],
				"' is for wrong run", NULL);
	       return TCL_ERROR;
	     }
	 }
       else
	 {
	   // If none was specified, use the current one.
	   (*fp[i])->currentAstrometricCalibration(astromH[i], oocRead);
	   astrom[i] = &astromH[i];
	 }

       // Get handle to the final photometric calibration
       sprintf(buffer, "-photo%d", i+1);
       if (ftcl_ArgIsPresent(argc, argv, buffer, runMergeQATbl) == 1)
	 {
	   // If one was specified, use it.
	   if (shTclAddrGetFromName (interp, finalName[i], (void **)&final[i],
				     "ooHandle(arFinalPhotoCalib)")
	       != TCL_OK)
	     {
	       Tcl_AppendResult(interp, runMergeQACmd,
				": error getting handle to final photometric calibration run '",
				finalName[i], "'", NULL);
	       return TCL_ERROR;
	     }
	   if ((*final[i])->framesPipelineRun() != (*fp[i]))
	     {
	       Tcl_AppendResult(interp, runMergeQACmd,
				": photometric calibration '",finalName[i],
				"' is for wrong run", NULL);
	       return TCL_ERROR;
	     }
	 }
       else
	 {
	   // If none was specified, use the current one.
	   (*fp[i])->currentPhotometricCalibration(finalH[i], oocRead);
	   if (finalH[i].isValid() == oocFalse)
	     {
	       Tcl_AppendResult(interp, runMergeQACmd, ": run '", fpName[i],
				"' has not been through final photometric calibration", NULL);
	       return TCL_ERROR;
	     }
	   final[i] = &finalH[i];
	 }

       // Get the exposure time (seconds)
       // HARDWIRE: All photometric chips have 2048 physical rows.
       int nRowsPhysical = 2048;
       expTime[i] = nRowsPhysical * (*fp[i])->imagingRun()->c_obs() * 1.e-6;
     }

   // Great circles must be for the same equinox.
   if (match)
     if (abs((*astrom[0])->equinox() - (*astrom[1])->equinox()) > 0.000001)
       {
	 Tcl_AppendResult(interp, runMergeQACmd, ": mismatched equinoxes",
			  NULL);
	 return TCL_ERROR;
       }

   // Set up cosmic color for use by atTransApply().
   taTransCosmicSet(cosmicColor, cosmicScatter);

   // Get/verify field range
   int field0 = (*fp[0])->field0();
   int field1 = field0 + (*fp[0])->nFields() - 1;
   if (firstField > -1)
     {
       if (firstField < field0 || firstField > field1)
	 {
	   Tcl_AppendResult(interp, selfMergeQACmd,
			    ": 'firstField' outside range of fields for this run",
			    NULL);
	   return TCL_ERROR;
	 }
     }
   else
     {
       firstField = field0;
     }
   if (lastField > -1)
     {
       if (lastField < field0 || lastField > field1)
	 {
	   Tcl_AppendResult(interp, selfMergeQACmd,
			    ": 'lastField' outside range of fields for this run",
			    NULL);
	   return TCL_ERROR;
	 }
     }
   else
     {
       lastField = field1;
     }
   int nFields = lastField - firstField + 1;

   // Fetch the database containing the objects for the second run, and the
   // camera column.
   ooHandle(ooDBObj) db2;
   uint8 camCol = 0;
   if (match)
     {
       (*fp[1])->objectList((*fp[1])->field0())->objectCont().containedIn(db2);
       camCol = (*fp[1])->camCol();
     }

   // Allocate a CHAIN to return the matching statistics for each field in the
   // first run
   CHAIN *ch = shChainNewBlock("TSMATCHSTATS", nFields);

   // If object chains are desired, allocate them.
   CHAIN *starChain = NULL;
   CHAIN *galaxyChain = NULL;
   CHAIN *otherChain = NULL;
   CHAIN *fieldChain1 = NULL;
   CHAIN *fieldChain2 = NULL;
   CHAIN *missingChain = NULL;
   if (objChains)
     {
       if (match)
	 {
	   starChain = shChainNew("TSMATCHEDPAIR");
	   galaxyChain = shChainNew("TSMATCHEDPAIR");
	   otherChain = shChainNew("TSMATCHEDPAIR");
	 }
       if (full)
	 {
	   fieldChain1 = shChainNew("TA_FIELD_INFO");
	   if (match)
	     {
	       fieldChain2 = shChainNew("TA_FIELD_INFO");
	     }
	   else
	     {
	       starChain = shChainNew("TA_CALIB_OBJ");
	       galaxyChain = shChainNew("TA_CALIB_OBJ");
	       otherChain = shChainNew("TA_CALIB_OBJ");
	     }
	 }
     }
   if (missingDistance >= 0) missingChain = shChainNew("TA_CALIB_OBJ");

   // Create field info chain and array for second run
   TA_FIELD_INFO *otherFieldArray[1500];
   if (full && match)
     {
       // Set up array of camRow versus filter, in same order as filter order
       // in array fields in TA_CALIB_OBJ (ugriz).  This is used when getting
       // TRANS structures for specific filters.
       const int camRow[5] = {3, 5, 1, 2, 4};

       int otherField0 = (*fp[1])->field0();
       int otherField1 = otherField0 + (*fp[1])->nFields() - 1;
       for (int32 field = otherField0; field <= otherField1; field++)
	 {
	   TA_FIELD_INFO *fieldInfo =
	     (TA_FIELD_INFO*) shMalloc(sizeof(TA_FIELD_INFO));
	   shAssert(fieldInfo != NULL);
	   if (shChainElementAddByPos(fieldChain2, fieldInfo, "TA_FIELD_INFO",
				      TAIL, AFTER)
	       != SH_SUCCESS)
	     {
	       Tcl_AppendResult(interp, runMergeQACmd,
				": couldn't append field info chain", NULL);
	       return TCL_ERROR;
	     }
	   otherFieldArray[field] = fieldInfo;
	   fieldInfo->field = field;
	   ooHandle(arPSCalib) ps =
	     (*fp[1])->postageStampCalibration(oocRead);
	   arFieldPhotoTrans psTrans = ps->trans(field);
	   arFieldPhotoTrans2 psTrans2 = ps->trans2(field);
	   fieldInfo->psp_status = psTrans2.psp_status();
	       arFieldPTrans ptrans = (*final[1])->ptrans(field);
	   for (int kk = 0; kk < 5; kk++)
	     {
	       TRANS trans =
		 astromH[1]->astrotoolsTrans(camRow[kk]*10+camCol, field);
	       fieldInfo->mjd[kk] = trans.mjd;
	       fieldInfo->seeing[kk] =
		 psTrans.filter((arFilter)kk).psf_width();
	       fieldInfo->status[kk] = psTrans2.status((arFilter)kk);
	       fieldInfo->psf_ap_correctionErr[kk] =
		 psTrans2.psf_ap_correctionErr((arFilter)kk);

	       // Must calibrate sky from DNs/pixel to mag/arcsec^2
	       arPTrans dbPTrans = ptrans.filter((arFilter)kk);
	       double scale = at_deg2Asec * sqrt((double)(trans.b * trans.f -
							  trans.e * trans.c));
	       double scale2 = scale * scale;
	       ooHandle(arObjectList) olist = (*fp[1])->objectList(field);
	       double mag, magErr;
	       int status;
	       taLCaliApply(kk,0,olist->frameStat((arFilter)kk).sky_frames_sub() /
			    scale2, 0.,
			    olist->frameStat((arFilter)adj[kk]).sky_frames_sub() / scale2,
			    0., expTime[1], expTime[1], trans.airmass,
			    dbPTrans.a().value(), dbPTrans.k().value(),
			    dbPTrans.b().value(), dbPTrans.c().value(),
			    sign[kk], cosmicSkyColor[cosmicIdx[kk]],
			    cosmicSkyScatter[cosmicIdx[kk]], &mag, &magErr,
			    &status);
	       fieldInfo->sky_frames_sub[kk] = mag;
	     }
	 }
     }

   // For each field ...
   TSMATCHSTATS *stats;
   CURSOR_T cursor = shChainCursorNew(ch);
   for (int32 field = firstField; field <= lastField; field++)
     {
       printf("%4d", field);
       fflush(stdout);
       // Get the stats structure for this field
       if ((stats = (TSMATCHSTATS*) shChainWalk(ch, cursor, NEXT)) == NULL)
	 {
	   char fieldString[10];
	   sprintf(fieldString, "%d", field);
	   Tcl_AppendResult(interp, runMergeQACmd,
			    ": can't get stats for field ", fieldString, NULL);
	   return TCL_ERROR;
	 }
       stats->field = field;

       // Calculate the matching statistics for this field.
       if (tsFieldMatchStatsByArcsec(match,
				     (*fp[0])->objectList(field), *astrom[0],
				     *final[0], expTime[0], db2, camCol,
				     *astrom[1], *final[1], expTime[1],
				     minMag, maxMag, filter[0], mags,
				     flags[0][0][0], flags[0][1][0],
				     flags[0][0][1], flags[0][1][1],
				     flags[1][0][0], flags[1][1][0],
				     flags[1][0][1], flags[1][1][1],
				     flagsFilter[0], stats, starChain,
				     galaxyChain, otherChain, fieldChain1,
				     otherFieldArray, full, raw,
				     missingDistance, missingChain)
	   != SH_SUCCESS)
	 {
	   char fieldString[10];
	   sprintf(fieldString, "%d", field);
	   Tcl_AppendResult(interp, runMergeQACmd,
			    ": error in field ", fieldString, NULL);
	   return TCL_ERROR;
	 }
     }
   printf("\n");
   fflush(stdout);

   // Return a handle to the matching statistics chain
   char vName[HANDLE_NAMELEN];
   if (shTclHandleNew(interp, vName, "CHAIN", (void *)ch) != TCL_OK)
     {
       Tcl_AppendResult(interp, runMergeQACmd,
			": error allocating chain TCL handle", NULL);
       return TCL_ERROR;
     }
   Tcl_AppendElement(interp, vName);
   if ((match && objChains) || (! match && full))
     {
       if (shTclHandleNew(interp, vName, "CHAIN", (void *)starChain) != TCL_OK)
	 {
	   Tcl_AppendResult(interp, runMergeQACmd,
			    ": error allocating chain TCL handle", NULL);
	   return TCL_ERROR;
	 }
       Tcl_AppendElement(interp, vName);
       if (shTclHandleNew(interp, vName, "CHAIN", (void *)galaxyChain) != TCL_OK)
	 {
	   Tcl_AppendResult(interp, runMergeQACmd,
			    ": error allocating chain TCL handle", NULL);
	   return TCL_ERROR;
	 }
       Tcl_AppendElement(interp, vName);
       if (shTclHandleNew(interp, vName, "CHAIN", (void *)otherChain) !=TCL_OK)
	 {
	   Tcl_AppendResult(interp, runMergeQACmd,
			    ": error allocating chain TCL handle", NULL);
	   return TCL_ERROR;
	 }
       Tcl_AppendElement(interp, vName);
       if (full)
	 {
	   if (shTclHandleNew(interp, vName, "CHAIN", (void *)fieldChain1)
	       != TCL_OK)
	     {
	       Tcl_AppendResult(interp, runMergeQACmd,
				": error allocating chain TCL handle", NULL);
	       return TCL_ERROR;
	     }
	   Tcl_AppendElement(interp, vName);
	   if (match)
	     {
	       if (shTclHandleNew(interp, vName, "CHAIN", (void *)fieldChain2)
		   != TCL_OK)
		 {
		   Tcl_AppendResult(interp, runMergeQACmd,
				    ": error allocating chain TCL handle",
				    NULL);
		   return TCL_ERROR;
		 }
	       Tcl_AppendElement(interp, vName);
	     }
	 }
     }
   if (missingDistance >= 0)
     {
       if (shTclHandleNew(interp, vName, "CHAIN", (void *)missingChain)
	   != TCL_OK)
	 {
	   Tcl_AppendResult(interp, runMergeQACmd,
			    ": error allocating chain TCL handle", NULL);
	   return TCL_ERROR;
	 }
       Tcl_AppendElement(interp, vName);
     }
   return TCL_OK;
}

/*============================================================================
**
**
** TCL VERB: runObjectsGet
**
**
**	Fetch all OK_RUN objects in a run and return them as a TA_CALIB_OBJ
**	chain.  The objects are not calibrated, except for the astrometric
**	parameters (which are calculated using cosmic colors).  Objects to
**	consider are
**      restricted to those with specified flags turned on and specified flags
**      turned off (separate sets of flags for the two different runs); these
**      flags can be either the objc_flags or the flags in the specified
**      filter.  Flags are specified as a TCL list of ascii values.  For
**      example, to fetch only objects which are EDGE objects and not
**      blended (neither BLEND nor CHILD flags set):
**      
**          runObjectsGet ... -flagsOn EDGE -flagsOff {BLEND CHILD}
**      

**	The objects can be further restrict by specifying a standard
**	objectivity predicate.  Objects can further be restricted to be at
**	least
**	an "nSigma" detection in the r' band (as measured by its counts_model).
**      

**	A 2-element TCL list is returned, where the first element is the object
**	chain, and the second element is a field info (TA_FIELD_INFO) chain.
**
**
**
**============================================================================
*/
static char runObjectsGetCmd[] = "runObjectsGet";
static int runObjectsGetFlg = FTCL_ARGV_NO_LEFTOVERS;
static ftclArgvInfo runObjectsGetTbl[] = {
   {NULL, FTCL_ARGV_HELP, NULL, NULL,
    "Fetch uncalibrated objects from a run\n"},
   {"", FTCL_ARGV_STRING, NULL, NULL,
    "Handle of Frames Pipeline run (ooHandle(arFramesPipeRun))"},
   {"-firstField", FTCL_ARGV_INT, NULL, NULL, "First field"},
   {"-lastField", FTCL_ARGV_INT, NULL, NULL, "Last field"},
   {"-flagsFilter", FTCL_ARGV_STRING, NULL, NULL, "Filter to fetch flags from (u, g, r, i, or z --- defaults to objc_flags"},
   {"-flagsOn", FTCL_ARGV_STRING, NULL, NULL,
    "Flags which must be on for fetched objects"},
   {"-flagsOff", FTCL_ARGV_STRING, NULL, NULL,
    "Flags which must be off for fetched objects"},
   {"-predicate", FTCL_ARGV_STRING, NULL, NULL,
    "Standard objectivity query predicate"},
   {"-nSigma", FTCL_ARGV_DOUBLE, NULL, NULL,
    "Require objects to be an 'nSigma' detection in r' band (as measured by its counts_model)"},
   {NULL, FTCL_ARGV_END, NULL, NULL, NULL}
};

static int
tclRunObjectsGet (ClientData clientData,
		  Tcl_Interp *interp,
		  int argc,
		  char **argv)
{
   ooHandle(arFramesPipeRun) *fp = NULL;
   char *fpName = NULL;
   int firstField = -1, lastField = -1;
   char *flagsFilter = " ";
   char *flagsString[2] = {NULL, NULL};
   char *predicate = NULL;
   double nSigma = -1;
   int status;

   // Parse the arguments
   runObjectsGetTbl[1].dst = &fpName;
   runObjectsGetTbl[2].dst = &firstField;
   runObjectsGetTbl[3].dst = &lastField;
   runObjectsGetTbl[4].dst = &flagsFilter;
   runObjectsGetTbl[5].dst = &(flagsString[0]);
   runObjectsGetTbl[6].dst = &(flagsString[1]);
   runObjectsGetTbl[7].dst = &predicate;
   runObjectsGetTbl[8].dst = &nSigma;
   if ((status = shTclParseArgv(interp, &argc, argv, runObjectsGetTbl,
				runObjectsGetFlg, runObjectsGetCmd))
       != FTCL_ARGV_SUCCESS)
     return status;

   // Decode flags strings
   GenericEnum* enum1 = GenericEnum::Find("AR_DETECTION_FLAG");
   shAssert(enum1 != NULL);
   GenericEnum* enum2 = GenericEnum::Find("AR_DETECTION_FLAG2");
   shAssert(enum2 != NULL);
   char *prefix1 = "AR_DFLAG_";
   char *prefix2 = "AR_DFLAG2_";
   int flags[2][2] = {0, 0, 0, 0}; // [enum1/2][on/off]
   for (int j = 0; j < 2; j++)
     {
       if (flagsString[j] == NULL) continue;
       char **listArgv;
       int nFlags = 0;
       if (Tcl_SplitList(interp, flagsString[j], &nFlags, &listArgv)
	   != TCL_OK)
	 {
	   Tcl_AppendResult(interp, runObjectsGetCmd,
			    ": couldn't parse flags string", NULL);
	   return TCL_ERROR;
	 }
       char buffer[100];
       for (int k = 0; k < nFlags; k++)
	 {
	   if (strlen(listArgv[k]) > 40)
	     {
	       Tcl_AppendResult(interp, runObjectsGetCmd,
				": illegal flag '", listArgv[k], "'", NULL);
	       return TCL_ERROR;
	     }
	   sprintf(buffer, "%s%s", prefix1, listArgv[k]);
	   if (enum1->validKey(buffer) == TRUE)
	     flags[0][j] |= enum1->asInteger(buffer);
	   else
	     {
	       sprintf(buffer, "%s%s", prefix2, listArgv[k]);
	       if (enum2->validKey(buffer) == TRUE)
		 flags[1][j] |= enum2->asInteger(buffer);
	       else
		 {
		   Tcl_AppendResult(interp, runObjectsGetCmd,
				    ": illegal flag '", listArgv[k], "'",
				    NULL);
		   return TCL_ERROR;
		 }
	     }
	 }
     }

   // Get handle to the frame pipeline run
   if (shTclAddrGetFromName (interp, fpName, (void **)&fp,
			     "ooHandle(arFramesPipeRun)") != TCL_OK)
     {
       Tcl_AppendResult(interp, runObjectsGetCmd,
			": error getting handle to run '", fpName,
			",", NULL);
       return TCL_ERROR;
     }
   
   // Verify that this run has been through setObjectStatus
   if ((*fp)->setObjectStatus().isValid() == oocFalse)
     {
       Tcl_AppendResult(interp, runObjectsGetCmd, ": run '", fpName,
			"' has not been through setObjectStatus", NULL);
       return TCL_ERROR;
     }
   
   // Get astrometric calibration
   ooHandle(arAstromCalib) astromCalib =
     (*fp)->currentAstrometricCalibration();

   // Set up cosmic color for use by atTransApply().
   taTransCosmicSet(cosmicColor, cosmicScatter);

   // Get/verify field range
   int field0 = (*fp)->field0();
   int field1 = field0 + (*fp)->nFields() - 1;
   if (firstField > -1)
     {
       if (firstField < field0 || firstField > field1)
	 {
	   Tcl_AppendResult(interp, runObjectsGetCmd,
			    ": 'firstField' outside range of fields for this run",
			    NULL);
	   return TCL_ERROR;
	 }
     }
   else
     {
       firstField = field0;
     }
   if (lastField > -1)
     {
       if (lastField < field0 || lastField > field1)
	 {
	   Tcl_AppendResult(interp, runObjectsGetCmd,
			    ": 'lastField' outside range of fields for this run",
			    NULL);
	   return TCL_ERROR;
	 }
     }
   else
     {
       lastField = field1;
     }

   // Allocate chains
   CHAIN *objChain = shChainNew("TA_CALIB_OBJ");
   CHAIN *fieldChain = shChainNew("TA_FIELD_INFO");

   // For each field ...
   for (int32 field = firstField; field <= lastField; field++)
     {
       printf("%4d", field);
       fflush(stdout);

       // Calculate the matching statistics for this field.
       if (tsFieldObjectsGet((*fp)->objectList(field), astromCalib, nSigma,
			     flags[0][0], flags[1][0],
			     flags[0][1], flags[1][1],
			     flagsFilter[0], predicate,
			     objChain, fieldChain)
	   != SH_SUCCESS)
	 {
	   char fieldString[10];
	   sprintf(fieldString, "%d", field);
	   Tcl_AppendResult(interp, runObjectsGetCmd,
			    ": error in field ", fieldString, NULL);
	   return TCL_ERROR;
	 }
     }
   printf("\n");
   fflush(stdout);

   // Return a handle to the matching statistics chain
   char vName[HANDLE_NAMELEN];
   if (shTclHandleNew(interp, vName, "CHAIN", (void *)objChain) != TCL_OK)
     {
       Tcl_AppendResult(interp, runObjectsGetCmd,
			": error allocating chain TCL handle", NULL);
       return TCL_ERROR;
     }
   Tcl_AppendElement(interp, vName);
   if (shTclHandleNew(interp, vName, "CHAIN", (void *)fieldChain) != TCL_OK)
     {
       Tcl_AppendResult(interp, runObjectsGetCmd,
			": error allocating chain TCL handle", NULL);
       return TCL_ERROR;
     }
   Tcl_AppendElement(interp, vName);
   return TCL_OK;
}

/*============================================================================
**
**
** TCL VERB: matchedPairChainDel
**
**
**	Delete a TSMATCHPAIR chain and all its entries.
**
**
**
**============================================================================
*/
static char matchedPairChainDelCmd[] = "matchedPairChainDel";
static int matchedPairChainDelFlg = FTCL_ARGV_NO_LEFTOVERS;
static ftclArgvInfo matchedPairChainDelTbl[] = {
   {NULL, FTCL_ARGV_HELP, NULL, NULL,
    "Delete a TSMATCHPAIR chain and all its entries\n"},
   {"", FTCL_ARGV_STRING, NULL, NULL,
    "Handle to the chain"},
   {NULL, FTCL_ARGV_END, NULL, NULL, NULL}
};

static int
tclMatchedPairChainDel (ClientData clientData,
			Tcl_Interp *interp,
			int argc,
			char **argv)
{
  CHAIN *chain = NULL;
  char *chainName = NULL;
  CURSOR_T crsr;
  TSMATCHEDPAIR *elem = NULL;
  int status;

  // Parse the arguments
  matchedPairChainDelTbl[1].dst = &chainName;
  if ((status = shTclParseArgv(interp, &argc, argv, matchedPairChainDelTbl,
			       matchedPairChainDelFlg, matchedPairChainDelCmd))
      != FTCL_ARGV_SUCCESS)
    return status;

   
  // Translate the handle to the CHAIN
  if (shTclAddrGetFromName (interp, chainName, (void **)&chain,
			    "CHAIN") != TCL_OK)
    {
      Tcl_SetResult(interp,"arg does not appear to be a CHAIN!",TCL_STATIC);
      return(TCL_ERROR);
    }


  // Check that we've got the right CHAIN type
  if (shChainTypeGet(chain) != shTypeGetFromName ("TSMATCHEDPAIR"))
    {
      Tcl_SetResult(interp,"wrong chain type for chain!",TCL_STATIC);
      return(TCL_ERROR);
    }

  // Done with all the handles, etc. Now the action begins
  crsr = shChainCursorNew(chain);
  while ((elem = (TSMATCHEDPAIR*) shChainWalk(chain, crsr, NEXT)) != NULL)
    {
      if (elem->obj1 != NULL) shFree(elem->obj1);
      if (elem->obj2 != NULL) shFree(elem->obj2);
      shFree(elem);
    }
  shChainCursorDel(chain,crsr);
  shChainDel(chain);
  p_shTclHandleDel(interp, chainName); 

  Tcl_SetResult(interp,"",TCL_STATIC);
  return TCL_OK;
}

/*============================================================================
**
**
** TCL VERB: mergeQABin
**
**
**	Bin a TSMATCHSTATS chain by the specified number of entries.
**
**
**
**============================================================================
*/
static char mergeQABinCmd[] = "mergeQABin";
static int mergeQABinFlg = FTCL_ARGV_NO_LEFTOVERS;
static ftclArgvInfo mergeQABinTbl[] = {
   {NULL, FTCL_ARGV_HELP, NULL, NULL,
    "Bin a TSMATCHSTATS chain by the specified number of entries\n"},
   {"", FTCL_ARGV_STRING, NULL, NULL,
    "Handle to the TSMATCHSTATS chain"},
   {"", FTCL_ARGV_INT, NULL, NULL,
    "number of entries per bin"},
   {NULL, FTCL_ARGV_END, NULL, NULL, NULL}
};

static int
tclMergeQABin (ClientData clientData,
			Tcl_Interp *interp,
			int argc,
			char **argv)
{
  CHAIN *inChain = NULL;
  char *inChainName = NULL;
  int binSize = 0;
  int status;

  // Parse the arguments
  mergeQABinTbl[1].dst = &inChainName;
  mergeQABinTbl[2].dst = &binSize;
  if ((status = shTclParseArgv(interp, &argc, argv, mergeQABinTbl,
			       mergeQABinFlg, mergeQABinCmd))
      != FTCL_ARGV_SUCCESS)
    return status;

  // Translate the handle to the CHAIN
  if (shTclAddrGetFromName (interp, inChainName, (void **)&inChain,
			    "CHAIN") != TCL_OK)
    {
      Tcl_AppendResult(interp, runMergeQACmd,
		       ": arg 'chain' not a CHAIN", NULL);
      return(TCL_ERROR);
    }

  // Check that we've got the right CHAIN type and a valid bin size
  if (shChainTypeGet(inChain) != shTypeGetFromName ("TSMATCHSTATS"))
    {
      Tcl_AppendResult(interp, runMergeQACmd,
		       ": arg 'chain' wrong CHAIN type", NULL);
      return(TCL_ERROR);
    }
  if (binSize < 2)
    {
      Tcl_AppendResult(interp, runMergeQACmd,
		       ": arg 'binSize' must be > 1", NULL);
      return(TCL_ERROR);
    }

  // Allocate output chain
  CHAIN *outChain = shChainNew("TSMATCHSTATS");

  // Loop through the input chain, binning as we go.
  int j = 0;
  int nElems = shChainSize(inChain);
  TSMATCHSTATS *bin = NULL;
  TSMATCHSTATS *elem = NULL;
  CURSOR_T crsr = shChainCursorNew(inChain);
  while ((elem = (TSMATCHSTATS*) shChainWalk(inChain, crsr, NEXT)) != NULL)
    {
      if (j % binSize == 0)
	{
	  // If new bin, allocate new structure and initialize all to zero.
	  // If not enough elements left-over for a new bin, then break out.
	  if (nElems - j < binSize) break;
	  bin = (TSMATCHSTATS*) shMalloc(sizeof(TSMATCHSTATS));
	  bin->field = 0;
	  tsMatchStatsInit(bin);
	}

      // Update counts.
      bin->field += elem->field;
      bin->nOverlap += elem->nOverlap;
      bin->nOverlapBright += elem->nOverlapBright;
      bin->nMismatches += elem->nMismatches;
      bin->nMismatchesBright += elem->nMismatchesBright;
      bin->nUnknown += elem->nMismatches;
      bin->nUnknownBright += elem->nMismatchesBright;
      bin->stars.nObjects += elem->stars.nObjects;
      bin->stars.nObjectsBright += elem->stars.nObjectsBright;
      bin->stars.nMatches += elem->stars.nMatches;
      bin->stars.nMatchesBright += elem->stars.nMatchesBright;
      bin->stars.row.nPairs += elem->stars.row.nPairs;
      bin->stars.row.av += elem->stars.row.av * elem->stars.row.nPairs;
      bin->stars.row.rms += (elem->stars.row.rms * elem->stars.row.rms +
			     elem->stars.row.av * elem->stars.row.av) *
	elem->stars.row.nPairs;
      bin->stars.col.nPairs += elem->stars.col.nPairs;
      bin->stars.col.av += elem->stars.col.av * elem->stars.col.nPairs;
      bin->stars.col.rms += (elem->stars.col.rms * elem->stars.col.rms +
			     elem->stars.col.av * elem->stars.col.av) *
	elem->stars.col.nPairs;
      for (int i = 0; i < nFilters; i++)
	{
	  bin->stars.mag[i].nPairs += elem->stars.mag[i].nPairs;
	  bin->stars.mag[i].av +=
	    elem->stars.mag[i].av * elem->stars.mag[i].nPairs;
	  bin->stars.mag[i].rms +=
	    (elem->stars.mag[i].rms * elem->stars.mag[i].rms +
	     elem->stars.mag[i].av * elem->stars.mag[i].av) *
	    elem->stars.mag[i].nPairs;
	}
      for (i = 0; i < nColors; i++)
	{
	  bin->stars.color[i].nPairs += elem->stars.color[i].nPairs;
	  bin->stars.color[i].av +=
	    elem->stars.color[i].av * elem->stars.color[i].nPairs;
	  bin->stars.color[i].rms +=
	    (elem->stars.color[i].rms * elem->stars.color[i].rms +
	     elem->stars.color[i].av * elem->stars.color[i].av) *
	    elem->stars.color[i].nPairs;
	}
      bin->galaxies.nObjects += elem->galaxies.nObjects;
      bin->galaxies.nObjectsBright += elem->galaxies.nObjectsBright;
      bin->galaxies.nMatches += elem->galaxies.nMatches;
      bin->galaxies.nMatchesBright += elem->galaxies.nMatchesBright;
      bin->galaxies.row.nPairs += elem->galaxies.row.nPairs;
      bin->galaxies.row.av +=
	elem->galaxies.row.av * elem->galaxies.row.nPairs;
      bin->galaxies.row.rms +=
	(elem->galaxies.row.rms * elem->galaxies.row.rms +
	 elem->galaxies.row.av * elem->galaxies.row.av) *
	elem->galaxies.row.nPairs;
      bin->galaxies.col.nPairs += elem->galaxies.col.nPairs;
      bin->galaxies.col.av +=
	elem->galaxies.col.av * elem->galaxies.col.nPairs;
      bin->galaxies.col.rms +=
	(elem->galaxies.col.rms * elem->galaxies.col.rms +
	 elem->galaxies.col.av * elem->galaxies.col.av) *
	elem->galaxies.col.nPairs;
      for (i = 0; i < nFilters; i++)
	{
	  bin->galaxies.mag[i].nPairs += elem->galaxies.mag[i].nPairs;
	  bin->galaxies.mag[i].av +=
	    elem->galaxies.mag[i].av * elem->galaxies.mag[i].nPairs;
	  bin->galaxies.mag[i].rms +=
	    (elem->galaxies.mag[i].rms * elem->galaxies.mag[i].rms +
	     elem->galaxies.mag[i].av * elem->galaxies.mag[i].av) *
	    elem->galaxies.mag[i].nPairs;
	}
      for (i = 0; i < nColors; i++)
	{
	  bin->galaxies.color[i].nPairs += elem->galaxies.color[i].nPairs;
	  bin->galaxies.color[i].av +=
	    elem->galaxies.color[i].av * elem->galaxies.color[i].nPairs;
	  bin->galaxies.color[i].rms +=
	    (elem->galaxies.color[i].rms * elem->galaxies.color[i].rms +
	     elem->galaxies.color[i].av * elem->galaxies.color[i].av) *
	    elem->galaxies.color[i].nPairs;
	}

      // If bin is finished, finish off stats
      j++;
      if (j % binSize == 0)
	{
	  bin->field /= binSize;

	  // Final star statistics
	  if (bin->stars.row.nPairs > 0)
	    {
	      bin->stars.row.av /= bin->stars.row.nPairs;
	      double meanSq = bin->stars.row.rms / bin->stars.row.nPairs -
		bin->stars.row.av * bin->stars.row.av;
	      bin->stars.row.rms = (meanSq > 0. ? sqrt(meanSq) : 0.);
	    }
	  if (bin->stars.col.nPairs > 0)
	    {
	      bin->stars.col.av /= bin->stars.col.nPairs;
	      double meanSq = bin->stars.col.rms / bin->stars.col.nPairs -
		bin->stars.col.av * bin->stars.col.av;
	      bin->stars.col.rms = (meanSq > 0. ? sqrt(meanSq) : 0.);
	    }
	  for (i = 0; i < nFilters; i++)
	    if (bin->stars.mag[i].nPairs > 0)
	      {
		bin->stars.mag[i].av /= bin->stars.mag[i].nPairs;
		double meanSq =
		  bin->stars.mag[i].rms / bin->stars.mag[i].nPairs -
		  bin->stars.mag[i].av * bin->stars.mag[i].av;
		bin->stars.mag[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.);
	      }
	  for (i = 0; i < nColors; i++)
	    if (bin->stars.color[i].nPairs > 0)
	      {
		bin->stars.color[i].av /= bin->stars.color[i].nPairs;
		double meanSq =
		  bin->stars.color[i].rms / bin->stars.color[i].nPairs -
		  bin->stars.color[i].av * bin->stars.color[i].av;
		bin->stars.color[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.);
	      }

	  // Final galaxy statistics
	  if (bin->galaxies.row.nPairs > 0)
	    {
	      bin->galaxies.row.av /= bin->galaxies.row.nPairs;
	      double meanSq =
		bin->galaxies.row.rms / bin->galaxies.row.nPairs -
		bin->galaxies.row.av * bin->galaxies.row.av;
	      bin->galaxies.row.rms = (meanSq > 0. ? sqrt(meanSq) : 0.);
	    }
	  if (bin->galaxies.col.nPairs > 0)
	    {
	      bin->galaxies.col.av /= bin->galaxies.col.nPairs;
	      double meanSq =
		bin->galaxies.col.rms / bin->galaxies.col.nPairs -
		bin->galaxies.col.av * bin->galaxies.col.av;
	      bin->galaxies.col.rms = (meanSq > 0. ? sqrt(meanSq) : 0.);
	    }
	  for (i = 0; i < nFilters; i++)
	    if (bin->galaxies.mag[i].nPairs > 0)
	      {
		bin->galaxies.mag[i].av /= bin->galaxies.mag[i].nPairs;
		double meanSq =
		  bin->galaxies.mag[i].rms / bin->galaxies.mag[i].nPairs -
		  bin->galaxies.mag[i].av * bin->galaxies.mag[i].av;
		bin->galaxies.mag[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.);
	      }
	  for (i = 0; i < nColors; i++)
	    if (bin->galaxies.color[i].nPairs > 0)
	      {
		bin->galaxies.color[i].av /= bin->galaxies.color[i].nPairs;
		double meanSq =
		  bin->galaxies.color[i].rms / bin->galaxies.color[i].nPairs -
		  bin->galaxies.color[i].av * bin->galaxies.color[i].av;
		bin->galaxies.color[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.);
	      }

	  // Add to chain
	  if (shChainElementAddByPos(outChain, bin, "TSMATCHSTATS", TAIL,
				     AFTER) != SH_SUCCESS)
	    {
	      shFree(bin);
	      taGenericChainDestroy(outChain);
	      Tcl_AppendResult(interp, runMergeQACmd,
			       ": couldn't append to output chain", NULL);
	      return TCL_ERROR;
	    }
	}
    }
  shChainCursorDel(inChain, crsr);

  // Return a handle to the binned statistics chain
  char vName[HANDLE_NAMELEN];
  if (shTclHandleNew(interp, vName, "CHAIN", (void *)outChain) != TCL_OK)
    {
      // Error.  Delete the chain and return an error.
      taGenericChainDestroy(outChain);
      Tcl_AppendResult(interp, runMergeQACmd,
		       ": error allocating output chain TCL handle", NULL);
      return TCL_ERROR;
    }
  Tcl_SetResult(interp, vName, TCL_VOLATILE);
  return TCL_OK;
}

//****************************************************************************/
//
// Declare my new tcl verbs to tcl
//
void
tsTclMergeDeclare(Tcl_Interp *interp) 
{
  char *tclHelpFacil = "ts";
  shTclDeclare(interp, setObjectStatusDontCallCmd,
	       (Tcl_CmdProc *)tclSetObjectStatusDontCall,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, setObjectStatusDontCallTbl,
			       setObjectStatusDontCallFlg,
			       setObjectStatusDontCallCmd),
	       shTclGetUsage(interp, setObjectStatusDontCallTbl,
			     setObjectStatusDontCallFlg,
			     setObjectStatusDontCallCmd));
  shTclDeclare(interp, resetObjectStatusCmd,
	       (Tcl_CmdProc *)tclresetObjectStatus,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, resetObjectStatusTbl,
			       resetObjectStatusFlg,
			       resetObjectStatusCmd),
	       shTclGetUsage(interp, resetObjectStatusTbl,
			     resetObjectStatusFlg,
			     resetObjectStatusCmd));
  shTclDeclare(interp, runUnsetObjectStatusCmd,
	       (Tcl_CmdProc *)tclRunUnsetObjectStatus,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, runUnsetObjectStatusTbl,
			       runUnsetObjectStatusFlg,
			       runUnsetObjectStatusCmd),
	       shTclGetUsage(interp, runUnsetObjectStatusTbl,
			     runUnsetObjectStatusFlg,
			     runUnsetObjectStatusCmd));
  shTclDeclare(interp, selfMergeQACmd,
	       (Tcl_CmdProc *)tclSelfMergeQA,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, selfMergeQATbl, selfMergeQAFlg,
			       selfMergeQACmd),
	       shTclGetUsage(interp, selfMergeQATbl, selfMergeQAFlg,
			     selfMergeQACmd));
  shTclDeclare(interp, fieldPlotCmd,
	       (Tcl_CmdProc *)tclFieldPlot,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, fieldPlotTbl, fieldPlotFlg,
			       fieldPlotCmd),
	       shTclGetUsage(interp, fieldPlotTbl, fieldPlotFlg,
			     fieldPlotCmd));
  shTclDeclare(interp, runMergeDontCallCmd,
	       (Tcl_CmdProc *)tclRunMergeDontCall,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, runMergeDontCallTbl,
			       runMergeDontCallFlg, runMergeDontCallCmd),
	       shTclGetUsage(interp, runMergeDontCallTbl, runMergeDontCallFlg,
			     runMergeDontCallCmd));
  shTclDeclare(interp, runUnmergeCmd,
	       (Tcl_CmdProc *)tclRunUnmerge,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, runUnmergeTbl, runUnmergeFlg,
			       runUnmergeCmd),
	       shTclGetUsage(interp, runUnmergeTbl, runUnmergeFlg,
			     runUnmergeCmd));
  shTclDeclare(interp, runMergeQACmd,
	       (Tcl_CmdProc *)tclRunMergeQA,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, runMergeQATbl, runMergeQAFlg,
			       runMergeQACmd),
	       shTclGetUsage(interp, runMergeQATbl, runMergeQAFlg,
			     runMergeQACmd));
  shTclDeclare(interp, runObjectsGetCmd,
	       (Tcl_CmdProc *)tclRunObjectsGet,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, runObjectsGetTbl, runObjectsGetFlg,
			       runObjectsGetCmd),
	       shTclGetUsage(interp, runObjectsGetTbl, runObjectsGetFlg,
			     runObjectsGetCmd));
  shTclDeclare(interp, matchedPairChainDelCmd,
	       (Tcl_CmdProc *)tclMatchedPairChainDel,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, matchedPairChainDelTbl,
			       matchedPairChainDelFlg,
			       matchedPairChainDelCmd),
	       shTclGetUsage(interp, matchedPairChainDelTbl,
			     matchedPairChainDelFlg,
			     matchedPairChainDelCmd));
  shTclDeclare(interp, mergeQABinCmd,
	       (Tcl_CmdProc *)tclMergeQABin,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, mergeQABinTbl, mergeQABinFlg,
			       mergeQABinCmd),
	       shTclGetUsage(interp, mergeQABinTbl, mergeQABinFlg,
			     mergeQABinCmd));
}