/***************************************************************************** ****************************************************************************** ** ** ** FILE: merge.cc ** ** ** This file contains functions to merge, unmerge, and generate ** matching statistics. ** ** ** ** ** ENVIRONMENT: ** C++. ** ** REQUIRED PRODUCTS: ** photoSch ** ****************************************************************************** ****************************************************************************** */ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "tsMatchStats.h" #include "tsMerge.h" #include extern "C" { #include } #include #include #include #include "tsTarget.h" // Forward declarations static int cmpRow (const void *d1, const void *d2); int tsObjArrayBinomialSearch(double row, TSOBJ *array, int nObjects); static int cmpMu (const void *d1, const void *d2); int tsCalibObjArrayBinomialSearch(double mu, TSCALIBOBJ *array, int nObjects); int tsCalibObjArrayCorrelatedSearch(double mu, TSCALIBOBJ *array, int nObjects, int guess); // Number of filters const int nFilters = 5; const int nColors = 4; // Luptitude const double luptitudeA = -2.5 / log(10.0); const double luptitudeB = log(2.0 / LUPTITUDE_SCALE); double lasinh(double x) { return log(x + sqrt(1 + x * x)); } double luptitude(double x) { return luptitudeA * (lasinh(x / LUPTITUDE_SCALE) - luptitudeB); } /*============================================================================ ** ** ** ROUTINE: tsFieldUnsetObjectStatus ** ** ** Unset all object status bits for all objects in a single ** field of a Frames Pipeline Run. This also removes all links for all ** objects in the field to other objects belonging to the run. It does ** not remove links to objects in other Frames Pipeline Runs. Returns ** SH_SUCCESS if successful, SH_GENERIC_ERROR if an error occurred. ** ** **============================================================================ */ RET_CODE tsFieldUnsetObjectStatus(ooHandle(arObjectList)& field /* MODIFIED: the object list for a field */) { // Initialize the iterator over objects ooItr(arObject) itr; if (field->objects(itr, oocUpdate) == oocError) return SH_GENERIC_ERROR; // For each object ... while (itr.next()) { // Object has duplicates only if DUPLICATE flag is set. If it is, then // remove the links to the duplicates. if (itr->status(0) & AR_OBJECT_STATUS_DUPLICATE) if (itr->delDuplicates() == oocError) return SH_GENERIC_ERROR; // Set the object status to undetermined. itr->setStatus(0, 0); itr->setStatus(1, 0); } // Done return SH_SUCCESS; } /*============================================================================ ** ** ** ROUTINE: tsFieldUnmerge ** ** ** Remove links from all objects in a single field of a Frames Pipeline ** Run to objects in other Frames Pipeline Runs. It does not remove links ** to objects in the same Frames Pipeline Run. Returns SH_SUCCESS if ** successful, SH_GENERIC_ERROR if an error occurs. ** ** ** **============================================================================ */ RET_CODE tsFieldUnmerge(ooHandle(arObjectList)& field /* MODIFIED: the object list for a field */) { // Initialize the iterator over objects ooItr(arObject) itr; if (field->objects(itr, oocUpdate) == oocError) return SH_GENERIC_ERROR; // For each object ... while (itr.next()) { // Only OK_RUN objects were matched, so skip if not OK_RUN. // Remove all matches. if (itr->status(0) & AR_OBJECT_STATUS_OK_RUN) if (itr->delMatches() == oocError) return SH_GENERIC_ERROR; } // Done return SH_SUCCESS; } /*============================================================================ ** ** ** ROUTINE: tsMatchStatsInit ** ** ** Initialize a TSMATCHSTATS structure, setting all fields to 0. ** ** ** **============================================================================ */ void tsMatchStatsInit(TSMATCHSTATS *stats /* MODIFIED: matching statistics */) { // Initialize the statistics stats->nOverlap = 0; stats->nOverlapBright = 0; stats->nMismatches = 0; stats->nMismatchesBright = 0; stats->nUnknown = 0; stats->nUnknownBright = 0; stats->stars.nObjects = 0; stats->stars.nObjectsBright = 0; stats->stars.nMatches = 0; stats->stars.nMatchesBright = 0; stats->stars.row.av = 0.; stats->stars.row.rms = 0.; stats->stars.row.nPairs = 0.; stats->stars.col.av = 0.; stats->stars.col.rms = 0.; stats->stars.col.nPairs = 0.; for (int i = 0; i < nFilters; i++) { stats->stars.mag[i].av = 0; stats->stars.mag[i].rms = 0; stats->stars.mag[i].nPairs = 0; } for (i = 0; i < nColors; i++) { stats->stars.color[i].av = 0; stats->stars.color[i].rms = 0; stats->stars.color[i].nPairs = 0; } stats->galaxies.nObjects = 0; stats->galaxies.nObjectsBright = 0; stats->galaxies.nMatches = 0; stats->galaxies.nMatchesBright = 0; stats->galaxies.row.av = 0.; stats->galaxies.row.rms = 0.; stats->galaxies.row.nPairs = 0.; stats->galaxies.col.av = 0.; stats->galaxies.col.rms = 0.; stats->galaxies.col.nPairs = 0.; for (i = 0; i < nFilters; i++) { stats->galaxies.mag[i].av = 0; stats->galaxies.mag[i].rms = 0; stats->galaxies.mag[i].nPairs = 0.; } for (i = 0; i < nColors; i++) { stats->galaxies.color[i].av = 0; stats->galaxies.color[i].rms = 0; stats->galaxies.color[i].nPairs = 0; } } /*============================================================================ ** ** ** ROUTINE: tsFieldMatchStatsByPixels ** ** ** Calculate matching statistics between two adjacent fields in the same ** run. Detailed statistics are returned ** in an argument. Objects in the first "field" are assumed to be ** located at lower row number, by exactly "rowsPerField" rows, compared ** to objects in the "other" field. This is typically used to compare ** adjacent fields within a run, where the field number of "field" is ** always one larger than the field number for "other". ** The returned statistics give the row and ** column statistics in pixels, and the uncalibrated PSF magnitude ** statistics as magnitudes. Only GOOD objects of class STAR, GALAXY, or ** UNK are counted. ** ** ** **============================================================================ */ RET_CODE tsFieldMatchStatsByPixel(const ooHandle(arObjectList)& field /* IN: the object list for a field */, const ooHandle(arObjectList)& other /* IN: object list for the field to match against */, int rowsPerField, /* number of non-overlap rows per field */ double nSigma /* IN: Compile matching statistics for * nSigma detections (PSF, r') */, TSMATCHSTATS *stats /* MODIFIED: matching statistics */, int noBlends /* IN: don't included blended object */) { // Allocate handles for an iterator for, and the containing container // for, matching objects. ooItr(arObject) matches; ooHandle(ooContObj) contH; ooHandle(arObjectCont) otherContH; otherContH = other->objectCont(otherContH); // Initialize the statistics tsMatchStatsInit(stats); // If both fields are the same, then there will be no matches (the // calling procedure does this to compile the non-matching stats for the // first field in the run). int noMatches = (field == other); // Initialize the iterator over objects ooItr(arObject) itr; if (field->objects(itr, oocReadOnly) == oocError) return SH_GENERIC_ERROR; // For each object ... while (itr.next()) { // Skip BAD objects int status = itr->status(0); if (! (status & AR_OBJECT_STATUS_GOOD)) continue; // Count only stars, galaxies, and unknown arObjectType objType = itr->objc_type(); TSSTATS *matchStats; arFloatErr counts; if (objType == AR_OBJECT_TYPE_STAR) { matchStats = &stats->stars; matchStats->nObjects++; counts = itr->psfCounts(AR_R); } else if (objType == AR_OBJECT_TYPE_GALAXY) { matchStats = &stats->galaxies; matchStats->nObjects++; counts = itr->petroCounts(AR_R); } else if (objType == AR_OBJECT_TYPE_UNK) { matchStats = NULL; stats->nUnknown++; counts = itr->psfCounts(AR_R); } else continue; // Is this a "bright" object? int bright = (counts.value() > nSigma * counts.err() && counts.err() > 0.); if (bright) { if (objType == AR_OBJECT_TYPE_UNK) stats->nUnknownBright++; else matchStats->nObjectsBright++; } // If not a DUPLICATE object then there are no matching objects so // skip the matching test. if (! (status & AR_OBJECT_STATUS_DUPLICATE)) { if (! (status & AR_OBJECT_STATUS_OK_RUN)) { stats->nOverlap++; if (bright) stats->nOverlapBright++; } continue; } // Skip remaining if unknown, as we don't compile matching stats for // type unknown. if (objType == AR_OBJECT_TYPE_UNK) continue; // If no matches possible skip to next object. if (noMatches) continue; // Skip object if blended and we're skipping blends if (noBlends) if (itr->objc_flags() & AR_DFLAG_BLENDED || itr->objc_flags() & AR_DFLAG_CHILD) continue; // Initialize iterator over all matching objects in adjacent fields. if (itr->duplicates(matches, oocNoOpen) == oocError) return SH_GENERIC_ERROR; // For each matching object ... while (matches.next()) { // Fetch the container for this object. matches.containedIn(contH); // If the matching object is in the other field, increment // the matching statistics. if (contH == otherContH) { matchStats->nMatches++; if (objType != matches->objc_type()) { stats->nMismatches++; if (bright) stats->nMismatchesBright++; } // Add to statistics only if object is "bright" (at least an // nSigma detection). if (bright) { matchStats->nMatchesBright++; double rowc = itr->objc_rowc().value() + rowsPerField - matches->objc_rowc().value(); double colc = itr->objc_colc().value() - matches->objc_colc().value(); matchStats->row.av += rowc; matchStats->row.rms += rowc * rowc; matchStats->row.nPairs++; matchStats->col.av += colc; matchStats->col.rms += colc * colc; matchStats->col.nPairs++; arFloatErr bModel1, bModel2; for (int i = 0; i < nFilters; i++) { arFloatErr counts1, counts2, model1, model2; if (objType == AR_OBJECT_TYPE_STAR) { // Use PSF magnitudes for stars counts1 = itr->psfCounts((arFilter)i); counts2 = matches->psfCounts((arFilter)i); } else { // Use petrosian magnitudes for galaxies counts1 = itr->petroCounts((arFilter)i); counts2 = matches->petroCounts((arFilter)i); } if (counts1.err() > 0 && counts2.err() > 0) { double magDiff = luptitude(counts1.value()) - luptitude(counts2.value()); matchStats->mag[i].av += magDiff; matchStats->mag[i].rms += magDiff * magDiff; matchStats->mag[i].nPairs++; } // Update color stats model1 = itr->counts_model((arFilter)i); model2 = matches->counts_model((arFilter)i); if (i > 0) if (model1.err() > 0 && model2.err() > 0 && bModel1.err() > 0 && bModel2.err() > 0) { double magDiff = (luptitude(bModel1.value()) - luptitude(model1.value())) - (luptitude(bModel2.value()) - luptitude(model2.value())); matchStats->color[i-1].av += magDiff; matchStats->color[i-1].rms += magDiff * magDiff; matchStats->color[i-1].nPairs++; } bModel1 = model1; bModel2 = model2; } } } } } // Final star statistics if (stats->stars.row.nPairs > 0) { stats->stars.row.av /= stats->stars.row.nPairs; double meanSq = stats->stars.row.rms / stats->stars.row.nPairs - stats->stars.row.av * stats->stars.row.av; stats->stars.row.rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } if (stats->stars.col.nPairs > 0) { stats->stars.col.av /= stats->stars.col.nPairs; double meanSq = stats->stars.col.rms / stats->stars.col.nPairs - stats->stars.col.av * stats->stars.col.av; stats->stars.col.rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } for (int i = 0; i < nFilters; i++) if (stats->stars.mag[i].nPairs > 0) { stats->stars.mag[i].av /= stats->stars.mag[i].nPairs; double meanSq = stats->stars.mag[i].rms / stats->stars.mag[i].nPairs - stats->stars.mag[i].av * stats->stars.mag[i].av; stats->stars.mag[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } for (i = 0; i < nColors; i++) if (stats->stars.color[i].nPairs > 0) { stats->stars.color[i].av /= stats->stars.color[i].nPairs; double meanSq = stats->stars.color[i].rms / stats->stars.color[i].nPairs - stats->stars.color[i].av * stats->stars.color[i].av; stats->stars.color[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } // Final galaxy statistics if (stats->galaxies.row.nPairs > 0) { stats->galaxies.row.av /= stats->galaxies.row.nPairs; double meanSq = stats->galaxies.row.rms / stats->galaxies.row.nPairs - stats->galaxies.row.av * stats->galaxies.row.av; stats->galaxies.row.rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } if (stats->galaxies.col.nPairs > 0) { stats->galaxies.col.av /= stats->galaxies.col.nPairs; double meanSq = stats->galaxies.col.rms / stats->galaxies.col.nPairs - stats->galaxies.col.av * stats->galaxies.col.av; stats->galaxies.col.rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } for (i = 0; i < nFilters; i++) if (stats->galaxies.mag[i].nPairs > 0) { stats->galaxies.mag[i].av /= stats->galaxies.mag[i].nPairs; double meanSq = stats->galaxies.mag[i].rms / stats->galaxies.mag[i].nPairs - stats->galaxies.mag[i].av * stats->galaxies.mag[i].av; stats->galaxies.mag[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } for (i = 0; i < nColors; i++) if (stats->galaxies.color[i].nPairs > 0) { stats->galaxies.color[i].av /= stats->galaxies.color[i].nPairs; double meanSq = stats->galaxies.color[i].rms / stats->galaxies.color[i].nPairs - stats->galaxies.color[i].av * stats->galaxies.color[i].av; stats->galaxies.color[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } // Return the number of good matches. return SH_SUCCESS; } /*============================================================================ ** ** ** ROUTINE: tsFieldMatchStatsByArcsec ** ** ** Calculate matching statistics between a field and another Frames ** Pipeline run (specified as a handle to the database containing its ** objects). The astrometric ** and photometric calibrations to use for both the field and the other ** run are specified as arguments. The row and column differences (in ** arcsecs) are in calibrated great circle coordinates, where the ** coordinate system of the input field 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. **

** The statistics are returned as a TSMATCHSTATS structure. If the ** "starChain", "galaxyChain", and "otherChain" pointers are non-NULL, ** then matched star-star and galaxy-galaxy pairs are returned on the ** the first two chains, and mismatched pairs and pairs ** containing UNK type on the last chain. If "missingChain" is non-NULL, ** 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 above behavior occurs if "match" is set to 1. If set to 0, then ** no matching occurs. Only the non-matching statistics are calculated, ** with the sample defined as described above. If "full" is 1, and ** "starChain", "galaxyChain", and "otherChain" pointers are non-NULL, ** then the calibrated objects (TA_CALIB_OBJ) are returned on those ** chains (of class STAR, GALAXY, and UNK respectively). **

** If "raw" is set to 1, 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. ** ** ** **============================================================================ */ RET_CODE tsFieldMatchStatsByArcsec(int match /* IN: 1=match against other run, 0=don't */, const ooHandle(arObjectList)& field /* IN: the object list for a field */, const ooHandle(arAstromCalib)& fieldCalib /* IN: the astrometric calibration to use for the field */, const ooHandle(arFinalPhotoCalib)& fieldFinal /* IN: the final photometric calibration to use to the field */, double fieldExpTime /* IN: exposure time for the field (seconds) */, const ooHandle(ooDBObj)& other /* IN: database with objects from other run */, uint8 otherCamCol /* IN: camera column for the other run */, const ooHandle(arAstromCalib)& otherCalib /* IN: the astrometric calibration to use for the other run */, const ooHandle(arFinalPhotoCalib)& otherFinal /* IN: the final photometric calibration to use to the other run */, double otherExpTime /* IN: exposure time for the other run (seconds) */, double minMag /* IN: minimum mag for "BRIGHT" objs*/, double maxMag /* IN: maximum mag for "BRIGHT" objs*/, char magFilter /* IN: filter for min/maxMag(ugriz)*/, int mags /* IN: type of magnitude to use when * selecting 'bright' objects * 0=psf for stars, petro for gals, * 1=psf mags for both, * 2=petrosian mags for both, * 3=model mags for both */, int flagsOn /* IN: require matched field * objects in 'flagChain' to * have these flags on */, int flags2On /* IN: require matched field * objects in 'flagChain' to * have these flags on */, int flagsOff /* IN: require matched field * objects in 'flagChain' to * have these flags off */, int flags2Off /* IN: require matched field * objects in 'flagChain' to * have these flags off */, int flagsOtherOn /* IN: require matched other * objects in 'flagChain' to * have these flags on */, int flags2OtherOn /* IN: require matched other * objects in 'flagChain' to * have these flags on */, int flagsOtherOff /* IN: require matched other * objects in 'flagChain' to * have these flags off */, int flags2OtherOff/* IN: require matched other * objects in 'flagChain' to * have these flags off */, char flagsFilter /* IN: filter to fetch flags from * u, g, r, i, z or ' ' for * object flags */, TSMATCHSTATS *stats /* MOD: matching statistics */, CHAIN *starChain /* MOD: matched stars chain (NULL to * ignore) */, CHAIN *galaxyChain /* MOD: matched galaxies chain * (NULL to ignore) */, CHAIN *otherChain /* MOD: mismatched and matched * unknowns chain * (NULL to ignore) */, CHAIN *fieldChain /* MOD: field info chain */, TA_FIELD_INFO **otherFieldArray /*IN:array of pointers to field info for other run*/, int full /* IN: 0=mag and position only for object * chains, 1=full objects */, int raw /* IN: 1=don't calibrate lengths, * magnitudes, surface brightesses, * and angles in full version of * objects; 0=do */, double missingDistance /* IN: (pixels) */, CHAIN *missingChain /* MOD: chain of unmatched bright objects */ ) { /*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}; // Get mag filter index int fil; switch (magFilter) { case 'u': fil = 0; break; case 'g': fil = 1; break; case 'r': fil = 2; break; case 'i': fil = 3; break; case 'z': fil = 4; break; default: shErrStackPush("tsFieldMatchStatsByArcsec: illegal magFilter"); return (SH_GENERIC_ERROR); } // Get flag filter index arFilter flagFil; switch (flagsFilter) { case 'u': flagFil = AR_U; break; case 'g': flagFil = AR_G; break; case 'r': flagFil = AR_R; break; case 'i': flagFil = AR_I; break; case 'z': flagFil = AR_Z; break; case ' ': flagFil = AR_B; break; default: shErrStackPush("tsFieldMatchStatsByArcsec: illegal flagsFilter"); return (SH_GENERIC_ERROR); } // Check that magnitude range specified correctly. if (minMag > maxMag) { double tmp = minMag; minMag = maxMag; maxMag = tmp; } // Allocate handles for an iterator for, and the containing container // for, matching objects. ooItr(arObject) matches; ooHandle(ooContObj) contH; ooHandle(ooDBObj) dbH; // Initialize the statistics tsMatchStatsInit(stats); // 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 refCamRow = 1; int camCol = field->framesPipelineRun()->camCol(); int ccd = 10 * refCamRow + camCol; int runId = field->framesPipelineRun()->run(); int rerun = field->framesPipelineRun()->rerun(); int otherCcd; if (match) otherCcd = 10 * refCamRow + otherCamCol; // 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}; // Fetch the calibrations for this field. TRANS trans = fieldCalib->astrotoolsTrans(ccd, field->field()); arFieldPTrans ptrans = fieldFinal->ptrans(field->field()); double node = fieldCalib->node(); double incl = fieldCalib->incl(); double otherNode, otherIncl; // Setup some stuff used for fetching missing objects // HARDWIRED: frame geometry int nCols = 2048; int nRows = 1361; TRANS missingTrans; int fetchMissing = 0; int missingLo = 0; int missingHi = 0; double beginMu, beginNu, endMu, endNu; if (match) { otherNode = otherCalib->node(); otherIncl = otherCalib->incl(); if (missingChain != NULL) { // Fetch TRANS for field in other run closest to this field // Get mu at middle of target field double mu, nu; atTransApply(&trans, 'r', nRows / 2., 0., nCols / 2., 0., NULL, NULL, &mu, NULL, &nu, NULL); // Convert to other run great circle coordinates (degrees) double ra, dec; atGCToEq(mu, nu, &ra, &dec, node, incl); atEqToGC(ra, dec, &mu, &nu, otherNode, otherIncl); // Determine closest field in other run for (int ifield = otherCalib->field0(); ifield <= otherCalib->field0() + otherCalib->nFields() - 1; ifield++) { TRANS otherTrans = otherCalib->astrotoolsTrans(otherCcd, ifield); atTransApply(&otherTrans, 'r', 0., 0., nCols / 2., 0., NULL, NULL, &beginMu, NULL, &beginNu, NULL); atTransApply(&otherTrans, 'r', nRows, 0., nCols / 2., 0., NULL, NULL, &endMu, NULL, &endNu, NULL); atBound(&beginMu, mu-180., mu+180.); atBound(&endMu, mu-180., mu+180.); if (mu >= beginMu && mu <= endMu) { fetchMissing = 1; missingTrans = otherTrans; if (ifield > otherCalib->field0()) missingLo = 1; if (ifield < otherCalib->field0()+otherCalib->nFields()-1) missingHi = 1; break; } } } } TRANS trans1[5]; TRANS *transPtr1[5]; TRANS photoTrans1[5]; TRANS *photoTransPtr1[5]; TA_PTRANS ptrans1; TA_FIELD_INFO *fieldInfo1 = NULL; if (full) { // Set up field info structure for output shAssert(fieldChain != NULL); fieldInfo1 = (TA_FIELD_INFO*) shMalloc(sizeof(TA_FIELD_INFO)); shAssert(fieldInfo1 != NULL); if (shChainElementAddByPos(fieldChain, fieldInfo1, "TA_FIELD_INFO", TAIL, AFTER) != SH_SUCCESS) { shErrStackPush("tsFieldMatchStatsByArcsec: couldn't append field info chain"); return SH_GENERIC_ERROR; } fieldInfo1->field = field->field(); ooHandle(arPSCalib) ps = field->framesPipelineRun()->postageStampCalibration(oocRead); arFieldPhotoTrans psTrans = ps->trans(field->field()); arFieldPhotoTrans2 psTrans2 = ps->trans2(field->field()); fieldInfo1->psp_status = psTrans2.psp_status(); // Get astrometric calibrations used when photo ran ooHandle(arAstromCalib) photoAstrom= field->framesPipelineRun()->astrometricCalibration(); shAssert(fabs(photoAstrom->node() - node) < 0.0001); shAssert(fabs(photoAstrom->incl() - incl) < 0.0001); // Set up astrometric transformations in each filter. for (int kk = 0; kk < 5; kk++) { trans1[kk] = fieldCalib->astrotoolsTrans(camRow[kk]*10+camCol, field->field()); transPtr1[kk] = &(trans1[kk]); photoTrans1[kk] = photoAstrom->astrotoolsTrans(camRow[kk]*10+camCol, field->field()); photoTransPtr1[kk] = &(photoTrans1[kk]); // Correct TRANS used by photo for rowOffset and // colOffset values determined by photo double rowOffset = field->frameStat((arFilter)kk).rowOffset(); double colOffset = field->frameStat((arFilter)kk).colOffset(); if (kk == 2) { // r' is always the reference filter, and its rowOffset and // colOffset should always be zero shAssert(fabs(rowOffset) < 0.00001); shAssert(fabs(colOffset) < 0.00001); } else { photoTrans1[kk].a += (-photoTrans1[kk].b * rowOffset - photoTrans1[kk].c * colOffset); photoTrans1[kk].d += (-photoTrans1[kk].e * rowOffset - photoTrans1[kk].f * colOffset); } // Get photometric transformations in each filter arPTrans dbPTrans1 = ptrans.filter((arFilter)kk); ptrans1.a[kk].val = dbPTrans1.a().value(); ptrans1.a[kk].err = dbPTrans1.a().err(); ptrans1.b[kk].val = dbPTrans1.b().value(); ptrans1.b[kk].err = dbPTrans1.b().err(); ptrans1.c[kk].val = dbPTrans1.c().value(); ptrans1.c[kk].err = dbPTrans1.c().err(); ptrans1.k[kk].val = dbPTrans1.k().value(); ptrans1.k[kk].err = dbPTrans1.k().err(); // Fill in filter part of field info structure fieldInfo1->seeing[kk] = psTrans.filter((arFilter)kk).psf_width(); fieldInfo1->mjd[kk] = trans1[kk].mjd; fieldInfo1->status[kk] = psTrans2.status((arFilter)kk); fieldInfo1->psf_ap_correctionErr[kk] = psTrans2.psf_ap_correctionErr((arFilter)kk); // Must calibrate sky from DNs/pixel to mag/arcsec^2 double scale = at_deg2Asec * sqrt((double)(trans1[kk].b * trans1[kk].f - trans1[kk].e * trans1[kk].c)); double scale2 = scale * scale; double mag, magErr; int status; taLCaliApply(kk,0,field->frameStat((arFilter)kk).sky_frames_sub() / scale2, 0., field->frameStat((arFilter)adj[kk]).sky_frames_sub() / scale2, 0., fieldExpTime, fieldExpTime, trans1[kk].airmass, ptrans1.a[kk].val, ptrans1.k[kk].val, ptrans1.b[kk].val, ptrans1.c[kk].val, sign[kk], cosmicSkyColor[cosmicIdx[kk]], cosmicSkyScatter[cosmicIdx[kk]], &mag, &magErr, &status); fieldInfo1->sky_frames_sub[kk] = mag; } } // Initialize the iterator over objects ooItr(arObject) itr; if (field->objects(itr, oocReadOnly) == oocError) return SH_GENERIC_ERROR; // For each object ... while (itr.next()) { // Skip BAD objects int status = itr->status(0); if (! (status & AR_OBJECT_STATUS_GOOD)) continue; // Skip objects not meeting the flagsON, flagsOFF criterion int flags, flags2; if (flagFil == AR_B) { flags = itr->objc_flags(); flags2 = itr->objc_flags2(); } else { flags = itr->data()->detection(flagFil).flags(); flags2 = itr->data()->detection(flagFil).flags2(); } if (((flags & flagsOn) != flagsOn) || ((flags2 & flags2On) != flags2On) || (flags & flagsOff) || (flags2 & flags2Off)) continue; // Count only OK_RUN stars, galaxies, and unknown (which are the only // types matched). arObjectType objType = itr->objc_type(); TSSTATS *matchStats; arFloatErr counts, adjCounts; if (objType == AR_OBJECT_TYPE_STAR) { matchStats = &stats->stars; if (status & AR_OBJECT_STATUS_OK_RUN) matchStats->nObjects++; // Use PSF magnitudes for stars if (mags == 0) { counts = itr->psfCounts((arFilter)fil); adjCounts = itr->psfCounts((arFilter)adj[fil]); } } else if (objType == AR_OBJECT_TYPE_GALAXY) { matchStats = &stats->galaxies; if (status & AR_OBJECT_STATUS_OK_RUN) matchStats->nObjects++; // Use petrosian magnitudes for galaxies if (mags == 0) { counts = itr->petroCounts((arFilter)fil); adjCounts = itr->petroCounts((arFilter)adj[fil]); } } else if (objType == AR_OBJECT_TYPE_UNK) { if (status & AR_OBJECT_STATUS_OK_RUN) stats->nUnknown++; // Use PSF magnitudes for unknown if (mags == 0) { counts = itr->psfCounts((arFilter)fil); adjCounts = itr->psfCounts((arFilter)adj[fil]); } } else continue; // If a specific magnitude is requested for use in defining a "bright" // object, regardless of type, get that magnitude. if (mags == 1) { counts = itr->psfCounts((arFilter)fil); adjCounts = itr->psfCounts((arFilter)adj[fil]); } else if (mags == 2) { counts = itr->petroCounts((arFilter)fil); adjCounts = itr->petroCounts((arFilter)adj[fil]); } else if (mags == 3) { counts = itr->counts_model((arFilter)fil); adjCounts = itr->counts_model((arFilter)adj[fil]); } // Is this a "bright" object? int bright = 0; if (counts.value() > 0 && counts.err() > 0) { double mag, magErr; int tStatus; arPTrans p = ptrans.filter((arFilter)fil); taLCaliApply(fil,0,counts.value(), counts.err(), adjCounts.value(), adjCounts.err(), fieldExpTime, fieldExpTime, trans.airmass, p.a().value(), p.k().value(), p.b().value(), p.c().value(), sign[fil], cosmicColor[cosmicIdx[fil]], cosmicScatter[cosmicIdx[fil]], &mag, &magErr, &tStatus); if (mag >= minMag && mag <= maxMag) bright = 1; } if (bright && (status & AR_OBJECT_STATUS_OK_RUN)) { if (objType == AR_OBJECT_TYPE_UNK) stats->nUnknownBright++; else matchStats->nObjectsBright++; } // Only OK_RUN objects were matched, so skip if this object isn't OK_RUN. if (! (status & AR_OBJECT_STATUS_OK_RUN)) { if (! (status & AR_OBJECT_STATUS_DUPLICATE)) { stats->nOverlap++; if (bright) stats->nOverlapBright++; } continue; } // Are we skipping matching? if (! match) { // Yes. If we're not putting out full objects, then we're done // with this object. If object isn't bright, then we're also // done with it. if (! full || ! bright) continue; // Calibrate object and add to appropriate chain CHAIN *pairChain = NULL; if (objType == AR_OBJECT_TYPE_STAR) pairChain = starChain; else if (objType == AR_OBJECT_TYPE_GALAXY) pairChain = galaxyChain; else pairChain = otherChain; if (pairChain == NULL) continue; TA_CALIB_OBJ *obj = (TA_CALIB_OBJ*) shMalloc(sizeof(TA_CALIB_OBJ)); if (shChainElementAddByPos(pairChain, obj, "TA_CALIB_OBJ", TAIL, AFTER) != SH_SUCCESS) { shErrStackPush("tsFieldMatchStatsByArcsec: couldn't append object chain"); shFree(obj); return SH_GENERIC_ERROR; } if (tsObjCalibrate(itr, obj, &ptrans1, (const TRANS **)transPtr1, (const TRANS **)photoTransPtr1, node, incl, 'r', fieldExpTime, raw, 0) != SH_SUCCESS) { shErrStackPush("tsFieldMatchStatsByArcsec: couldn't calibrate object"); return SH_GENERIC_ERROR; } obj->field = fieldInfo1; obj->run = runId; obj->camCol = camCol; obj->rerun = rerun; obj->fieldId = field->field(); // Skip to next object. continue; } // Initialize iterator over all matching objects. if (itr->matches(matches, oocNoOpen) == oocError) return SH_GENERIC_ERROR; int first = 1; // For each matching object ... 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; } double mu1; double nu1; // Fetch the container and database for this object. matches.containedIn(contH); contH.containedIn(dbH); // If the matching object is not in the other database, skip it. if (dbH != other) continue; // Skip objects not meeting the flagsON, flagsOFF criterion if (flagFil == AR_B) { flags = matches->objc_flags(); flags2 = matches->objc_flags2(); } else { flags = matches->data()->detection(flagFil).flags(); flags2 = matches->data()->detection(flagFil).flags2(); } if (((flags & flagsOtherOn) != flagsOtherOn) || ((flags2 & flags2OtherOn) != flags2OtherOn) || (flags & flagsOtherOff) || (flags2 & flags2OtherOff)) continue; // Update counts if (objType != matches->objc_type()) stats->nMismatches++; else if (objType != AR_OBJECT_TYPE_UNK) matchStats->nMatches++; // Add to statistics and object chains only if object is "bright". if (! bright) continue; // Update bright counts CHAIN *pairChain = NULL; if (objType != matches->objc_type()) { // Mismatch. Put on "other" chain. stats->nMismatchesBright++; pairChain = otherChain; } else if (objType == AR_OBJECT_TYPE_UNK) { // UNK. Put on "other" chain pairChain = otherChain; } else { // Matched stars or galaxies. Put on their chain. matchStats->nMatchesBright++; if (objType == AR_OBJECT_TYPE_STAR) pairChain = starChain; else pairChain = galaxyChain; } // Fetch calibrated great circle coordinates for target // object, if we haven't already. if (first) { first = 0; arSphericalCoord coord1 = itr->gcCoord(trans, &ptrans); mu1 = coord1.longitude(); nu1 = coord1.latitude(); atBound(&mu1, 0., 360.); } // Fetch calibrated great circle coordinates for matching object. ooHandle(arObjectList) otherFieldH = ((ooHandle(arObjectCont)&)contH)->objectList(); int otherField = otherFieldH->field(); ooHandle(arFramesPipeRun) otherFPR = otherFieldH->framesPipelineRun(); int otherRun = otherFPR->run(); int otherRerun = otherFPR->rerun(); TRANS otherTrans = otherCalib->astrotoolsTrans(otherCcd, otherField); arFieldPTrans otherPTrans = otherFinal->ptrans(otherField); arSphericalCoord coord2 = matches->gcCoord(otherTrans, &otherPTrans); double mu2 = coord2.longitude(); atBound(&mu2, 0., 360.); double nu2 = coord2.latitude(); // Convert other coordinates to target run great circle // coordinates. double ra, dec; atGCToEq(mu2, nu2, &ra, &dec, otherNode, otherIncl); atEqToGC(ra, dec, &mu2, &nu2, node, incl); atBound(&mu2, 0., 360.); // Increment position statistics double rowDiff = mu1 - mu2; if (rowDiff > 180.) rowDiff -= 360.; else if (rowDiff < -180) rowDiff += 360.; rowDiff *= at_deg2Asec; double colDiff = (nu1 - nu2) * at_deg2Asec; if (objType != AR_OBJECT_TYPE_UNK && objType == matches->objc_type()) { matchStats->row.av += rowDiff; matchStats->row.rms += rowDiff * rowDiff; matchStats->row.nPairs++; matchStats->col.av += colDiff; matchStats->col.rms += colDiff * colDiff; matchStats->col.nPairs++; } // Add individual pair to output pair chain TSMATCHEDPAIR *pair; if (pairChain != NULL) { pair = (TSMATCHEDPAIR*) shMalloc(sizeof(TSMATCHEDPAIR)); pair->objc_flags[0] = itr->objc_flags(); pair->objc_flags[1] = matches->objc_flags(); pair->objc_flags2[0] = itr->objc_flags2(); pair->objc_flags2[1] = matches->objc_flags2(); if (shChainElementAddByPos(pairChain, pair, "TSMATCHEDPAIR", TAIL, AFTER) != SH_SUCCESS) { shErrStackPush("tsFieldMatchStatsByArcsec: couldn't append object chain"); shFree(pair); return SH_GENERIC_ERROR; } if (full) { pair->obj1 = (TA_CALIB_OBJ*) shMalloc(sizeof(TA_CALIB_OBJ)); pair->obj2 = (TA_CALIB_OBJ*) shMalloc(sizeof(TA_CALIB_OBJ)); // Get astrometric calibration used when photo ran // in the other field ooHandle(arAstromCalib) otherPhotoAstrom = otherFPR->astrometricCalibration(); shAssert(fabs(otherPhotoAstrom->node() - otherNode) < 0.0001); shAssert(fabs(otherPhotoAstrom->incl() - otherIncl) < 0.0001); TRANS trans2[5]; TRANS *transPtr2[5]; TRANS photoTrans2[5]; TRANS *photoTransPtr2[5]; TA_PTRANS ptrans2; for (int kk = 0; kk < 5; kk++) { // Get astrometric calibrations in other field trans2[kk] = otherCalib->astrotoolsTrans(camRow[kk]*10 + otherCamCol, otherField); transPtr2[kk] = &(trans2[kk]); photoTrans2[kk] = otherPhotoAstrom->astrotoolsTrans(camRow[kk]*10 +otherCamCol, otherField); photoTransPtr2[kk] = &(photoTrans2[kk]); // Correct TRANS used by photo for rowOffset and // colOffset values determined by photo double rowOffset = otherFieldH->frameStat((arFilter)kk).rowOffset(); double colOffset = otherFieldH->frameStat((arFilter)kk).colOffset(); if (kk == 2) { // r' is always the reference filter, and // its rowOffset and // colOffset should always be zero shAssert(fabs(rowOffset) < 0.00001); shAssert(fabs(colOffset) < 0.00001); } else { photoTrans2[kk].a += (-photoTrans2[kk].b * rowOffset - photoTrans2[kk].c * colOffset); photoTrans2[kk].d += (-photoTrans2[kk].e * rowOffset - photoTrans2[kk].f * colOffset); } arPTrans dbPTrans2 = otherPTrans.filter((arFilter)kk); ptrans2.a[kk].val = dbPTrans2.a().value(); ptrans2.a[kk].err = dbPTrans2.a().err(); ptrans2.b[kk].val = dbPTrans2.b().value(); ptrans2.b[kk].err = dbPTrans2.b().err(); ptrans2.c[kk].val = dbPTrans2.c().value(); ptrans2.c[kk].err = dbPTrans2.c().err(); ptrans2.k[kk].val = dbPTrans2.k().value(); ptrans2.k[kk].err = dbPTrans2.k().err(); } if (tsObjCalibrate(itr, pair->obj1, &ptrans1, (const TRANS **)transPtr1, (const TRANS **)photoTransPtr1, node, incl, 'r', fieldExpTime, raw, 0) != SH_SUCCESS) { shErrStackPush("tsFieldMatchStatsByArcsec: couldn't calibrate object"); return SH_GENERIC_ERROR; } pair->obj1->field = fieldInfo1; pair->obj1->run = runId; pair->obj1->camCol = camCol; pair->obj1->rerun = rerun; pair->obj1->fieldId = field->field(); if (tsObjCalibrate(matches, pair->obj2, &ptrans2, (const TRANS **)transPtr2, (const TRANS **)photoTransPtr2, otherNode, otherIncl, 'r', otherExpTime, raw, 0) != SH_SUCCESS) { shErrStackPush("tsFieldMatchStatsByArcsec: couldn't calibrate object"); return SH_GENERIC_ERROR; } pair->obj2->field = otherFieldArray[otherField]; pair->obj2->run = otherRun; pair->obj2->camCol = otherCamCol; pair->obj2->rerun = otherRerun; pair->obj2->fieldId = otherField; } else { pair->obj1 = NULL; pair->obj2 = NULL; } pair->mu = mu1; pair->nu = nu1; pair->dmu = rowDiff; pair->dnu = colDiff; } // Increment mag statistics float bMag1, bMag1Err, bMag2, bMag2Err; for (int i = 0; i < nFilters; i++) { // Calibrate counts to magnitudes arFloatErr counts1, adj1, counts2, adj2; if (objType == AR_OBJECT_TYPE_STAR) { // Use PSF magnitudes for stars counts1 = itr->psfCounts((arFilter)i); adj1 = itr->psfCounts((arFilter)adj[i]); counts2 = matches->psfCounts((arFilter)i); adj2 = matches->psfCounts((arFilter)adj[i]); } else { // Use petrosian magnitudes for galaxies counts1 = itr->petroCounts((arFilter)i); adj1 = itr->petroCounts((arFilter)adj[i]); counts2 = matches->petroCounts((arFilter)i); adj2 = matches->petroCounts((arFilter)adj[i]); } // Update mag stats if all counts are good if (counts1.err() >= 0 && counts2.err() >= 0 && adj1.err() >= 0 && adj2.err() >= 0) { double mag1, mag1Err, mag2, mag2Err; int tStatus; arPTrans p1 = ptrans.filter((arFilter)i); arPTrans p2 = otherPTrans.filter((arFilter)i); taLCaliApply(i,0,counts1.value(), counts1.err(), adj1.value(), adj1.err(), fieldExpTime, fieldExpTime, trans.airmass, p1.a().value(), p1.k().value(), p1.b().value(), p1.c().value(), sign[i], cosmicColor[cosmicIdx[i]], cosmicScatter[cosmicIdx[i]], &mag1, &mag1Err, &tStatus); taLCaliApply(i,0,counts2.value(), counts2.err(), adj2.value(), adj2.err(), otherExpTime, otherExpTime, otherTrans.airmass, p2.a().value(), p2.k().value(), p2.b().value(), p2.c().value(), sign[i], cosmicColor[cosmicIdx[i]], cosmicScatter[cosmicIdx[i]], &mag2, &mag2Err, &tStatus); // Increment statistics double magDiff = mag1 - mag2; if (objType != AR_OBJECT_TYPE_UNK && objType == matches->objc_type()) { matchStats->mag[i].av += magDiff; matchStats->mag[i].rms += magDiff * magDiff; matchStats->mag[i].nPairs++; } if (pairChain != NULL) { pair->mag[i] = mag1; pair->dmag[i] = magDiff; } } // Get model mags for colors counts1 = itr->counts_model((arFilter)i); adj1 = itr->counts_model((arFilter)adj[i]); counts2 = matches->counts_model((arFilter)i); adj2 = matches->counts_model((arFilter)adj[i]); // If model counts are good, convert to mags if (counts1.err() >= 0 && counts2.err() >= 0 && adj1.err() >= 0 && adj2.err() >= 0) { double mag1, mag1Err, mag2, mag2Err; int tStatus; arPTrans p1 = ptrans.filter((arFilter)i); arPTrans p2 = otherPTrans.filter((arFilter)i); taLCaliApply(i,0,counts1.value(), counts1.err(), adj1.value(), adj1.err(), fieldExpTime, fieldExpTime, trans.airmass, p1.a().value(), p1.k().value(), p1.b().value(), p1.c().value(), sign[i], cosmicColor[cosmicIdx[i]], cosmicScatter[cosmicIdx[i]], &mag1, &mag1Err, &tStatus); taLCaliApply(i,0,counts2.value(), counts2.err(), adj2.value(), adj2.err(), otherExpTime, otherExpTime, otherTrans.airmass, p2.a().value(), p2.k().value(), p2.b().value(), p2.c().value(), sign[i], cosmicColor[cosmicIdx[i]], cosmicScatter[cosmicIdx[i]], &mag2, &mag2Err, &tStatus); // If possible, update color stats if (i > 0) if (bMag1Err >= 0 && bMag2Err >= 0) { // Increment color statistics double magDiff = (bMag1 - mag1) - (bMag2 - mag2); if (objType != AR_OBJECT_TYPE_UNK && objType == matches->objc_type()) { matchStats->color[i-1].av += magDiff; matchStats->color[i-1].rms += magDiff * magDiff; matchStats->color[i-1].nPairs++; } if (pairChain != NULL) { pair->color[i-1] = bMag1 - mag1; pair->dcolor[i-1] = magDiff; } } // Save blue mags for next band bMag1 = mag1; bMag1Err = mag1Err; bMag2 = mag2; bMag2Err = mag2Err; } else { // Undefined model mag. Save as undefined blue mags // for next band. bMag1Err = -1; bMag2Err = -1; } } } if (first == 1 && fetchMissing && bright) { // Didn't match anyone. Add to the missing chain. // Fetch calibrated great circle coordinates for target arSphericalCoord coord1 = itr->gcCoord(trans, &ptrans); double mu1 = coord1.longitude(); double nu1 = coord1.latitude(); // Convert to other run great circle coordinates double ra, dec; atGCToEq(mu1, nu1, &ra, &dec, node, incl); atEqToGC(ra, dec, &mu1, &nu1, otherNode, otherIncl); atBound(&mu1, beginMu - 180., beginMu + 180.); // Convert to pixel coordinates in other run nearest field double row, col; atTransInverseApply(&missingTrans, 'r', mu1, 0., nu1, 0., NULL, NULL, &row, NULL, &col, NULL); // Add to chain if within fair boundaries if (col > missingDistance && col < nCols - missingDistance) if ((row > missingDistance || missingLo) && (row < nRows - missingDistance || missingHi)) { TA_CALIB_OBJ *obj = (TA_CALIB_OBJ*) shMalloc(sizeof(TA_CALIB_OBJ)); if (shChainElementAddByPos(missingChain, obj, "TA_CALIB_OBJ", TAIL, AFTER) != SH_SUCCESS) { shErrStackPush("tsFieldMatchStatsByArcsec: couldn't append object chain"); shFree(obj); return SH_GENERIC_ERROR; } if (tsObjCalibrate(itr, obj, &ptrans1, (const TRANS **)transPtr1, (const TRANS **)photoTransPtr1, node, incl, 'r', fieldExpTime, raw, 0) != SH_SUCCESS) { shErrStackPush("tsFieldMatchStatsByArcsec: couldn't calibrate object"); return SH_GENERIC_ERROR; } obj->field = fieldInfo1; obj->run = runId; obj->camCol = camCol; obj->rerun = rerun; obj->fieldId = field->field(); } } } // Final star statistics if (stats->stars.row.nPairs > 0) { stats->stars.row.av /= stats->stars.row.nPairs; double meanSq = stats->stars.row.rms / stats->stars.row.nPairs - stats->stars.row.av * stats->stars.row.av; stats->stars.row.rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } if (stats->stars.col.nPairs > 0) { stats->stars.col.av /= stats->stars.col.nPairs; double meanSq = stats->stars.col.rms / stats->stars.col.nPairs - stats->stars.col.av * stats->stars.col.av; stats->stars.col.rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } for (int i = 0; i < nFilters; i++) if (stats->stars.mag[i].nPairs > 0) { stats->stars.mag[i].av /= stats->stars.mag[i].nPairs; double meanSq = stats->stars.mag[i].rms / stats->stars.mag[i].nPairs - stats->stars.mag[i].av * stats->stars.mag[i].av; stats->stars.mag[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } for (i = 0; i < nColors; i++) if (stats->stars.color[i].nPairs > 0) { stats->stars.color[i].av /= stats->stars.color[i].nPairs; double meanSq = stats->stars.color[i].rms / stats->stars.color[i].nPairs - stats->stars.color[i].av * stats->stars.color[i].av; stats->stars.color[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } // Final galaxy statistics if (stats->galaxies.row.nPairs > 0) { stats->galaxies.row.av /= stats->galaxies.row.nPairs; double meanSq = stats->galaxies.row.rms / stats->galaxies.row.nPairs - stats->galaxies.row.av * stats->galaxies.row.av; stats->galaxies.row.rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } if (stats->galaxies.col.nPairs > 0) { stats->galaxies.col.av /= stats->galaxies.col.nPairs; double meanSq = stats->galaxies.col.rms / stats->galaxies.col.nPairs - stats->galaxies.col.av * stats->galaxies.col.av; stats->galaxies.col.rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } for (i = 0; i < nFilters; i++) if (stats->galaxies.mag[i].nPairs > 0) { stats->galaxies.mag[i].av /= stats->galaxies.mag[i].nPairs; double meanSq = stats->galaxies.mag[i].rms / stats->galaxies.mag[i].nPairs - stats->galaxies.mag[i].av * stats->galaxies.mag[i].av; stats->galaxies.mag[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } for (i = 0; i < nColors; i++) if (stats->galaxies.color[i].nPairs > 0) { stats->galaxies.color[i].av /= stats->galaxies.color[i].nPairs; double meanSq = stats->galaxies.color[i].rms / stats->galaxies.color[i].nPairs - stats->galaxies.color[i].av * stats->galaxies.color[i].av; stats->galaxies.color[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } // Return the number of matches. return SH_SUCCESS; } /*============================================================================ ** ** ** ROUTINE: tsFieldObjectsGet ** ** ** Fetch all OK_RUN objects in a field and append them to 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. The objects can be further restrict by specifying a standard ** objectivity predicate. Objects may also be restricted to be at least ** an "nSigma" detection in the r' band (as measured by its counts_model). **

** A TA_FIELD_INFO structure is created for the field, appended to the ** field chain, and all returned objects in the field have a pointer set ** to it. ** ** ** **============================================================================ */ RET_CODE tsFieldObjectsGet(const ooHandle(arObjectList)& field /* IN: the object list for a field */, const ooHandle(arAstromCalib)& fieldCalib /* IN: the astrometric calibration to use for the field */, double nSigma /* IN: Get objects whose r' model counts is this nSigma of its error */, int flagsOn /* IN: require objects to * have these flags on */, int flags2On /* IN: require objects to * have these flags on */, int flagsOff /* IN: require objects to * have these flags off */, int flags2Off /* IN: require objects to * have these flags off */, char flagsFilter /* IN: filter to fetch flags from * u, g, r, i, z or ' ' for * object flags */, const char *predicate /* IN: query predicate */, CHAIN *objChain /* MOD: objects chain */, CHAIN *fieldChain /* MOD: field info chain */ ) { // Get flag filter index arFilter flagFil; switch (flagsFilter) { case 'u': flagFil = AR_U; break; case 'g': flagFil = AR_G; break; case 'r': flagFil = AR_R; break; case 'i': flagFil = AR_I; break; case 'z': flagFil = AR_Z; break; case ' ': flagFil = AR_B; break; default: shErrStackPush("tsFieldObjectsGet: illegal flagsFilter"); return (SH_GENERIC_ERROR); } int camCol = field->framesPipelineRun()->camCol(); int runId = field->framesPipelineRun()->run(); int rerun = field->framesPipelineRun()->rerun(); // 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}; // Fetch the astrometric calibrations for this field. double node = fieldCalib->node(); double incl = fieldCalib->incl(); TRANS trans1[5]; TRANS *transPtr1[5]; TRANS photoTrans1[5]; TRANS *photoTransPtr1[5]; TA_FIELD_INFO *fieldInfo1 = NULL; // Set up field info structure for output shAssert(fieldChain != NULL); fieldInfo1 = (TA_FIELD_INFO*) shMalloc(sizeof(TA_FIELD_INFO)); shAssert(fieldInfo1 != NULL); if (shChainElementAddByPos(fieldChain, fieldInfo1, "TA_FIELD_INFO", TAIL, AFTER) != SH_SUCCESS) { shErrStackPush("tsFieldObjectsGet: couldn't append field info chain"); return SH_GENERIC_ERROR; } fieldInfo1->field = field->field(); ooHandle(arPSCalib) ps = field->framesPipelineRun()->postageStampCalibration(oocRead); arFieldPhotoTrans psTrans = ps->trans(field->field()); arFieldPhotoTrans2 psTrans2 = ps->trans2(field->field()); fieldInfo1->psp_status = psTrans2.psp_status(); // Get astrometric calibrations used when photo ran ooHandle(arAstromCalib) photoAstrom= field->framesPipelineRun()->astrometricCalibration(); shAssert(fabs(photoAstrom->node() - node) < 0.0001); shAssert(fabs(photoAstrom->incl() - incl) < 0.0001); // Set up astrometric transformations in each filter. for (int kk = 0; kk < 5; kk++) { trans1[kk] = fieldCalib->astrotoolsTrans(camRow[kk]*10+camCol, field->field()); transPtr1[kk] = &(trans1[kk]); photoTrans1[kk] = photoAstrom->astrotoolsTrans(camRow[kk]*10+camCol, field->field()); photoTransPtr1[kk] = &(photoTrans1[kk]); // Correct TRANS used by photo for rowOffset and // colOffset values determined by photo double rowOffset = field->frameStat((arFilter)kk).rowOffset(); double colOffset = field->frameStat((arFilter)kk).colOffset(); if (kk == 2) { // r' is always the reference filter, and its rowOffset and // colOffset should always be zero shAssert(fabs(rowOffset) < 0.00001); shAssert(fabs(colOffset) < 0.00001); } else { photoTrans1[kk].a += (-photoTrans1[kk].b * rowOffset - photoTrans1[kk].c * colOffset); photoTrans1[kk].d += (-photoTrans1[kk].e * rowOffset - photoTrans1[kk].f * colOffset); } // Fill in filter part of field info structure fieldInfo1->seeing[kk] = psTrans.filter((arFilter)kk).psf_width(); fieldInfo1->mjd[kk] = trans1[kk].mjd; fieldInfo1->status[kk] = psTrans2.status((arFilter)kk); fieldInfo1->psf_ap_correctionErr[kk] = psTrans2.psf_ap_correctionErr((arFilter)kk); } // Initialize the iterator over objects ooItr(arObject) itr; if (field->objects(itr, oocReadOnly, oocAll, predicate) == oocError) return SH_GENERIC_ERROR; // For each object ... while (itr.next()) { // Use OK_RUN objects only int status = itr->status(0); if (! (status & AR_OBJECT_STATUS_OK_RUN)) continue; // Skip objects not meeting the flagsON, flagsOFF criterion int flags, flags2; if (flagFil == AR_B) { flags = itr->objc_flags(); flags2 = itr->objc_flags2(); } else { flags = itr->data()->detection(flagFil).flags(); flags2 = itr->data()->detection(flagFil).flags2(); } if (((flags & flagsOn) != flagsOn) || ((flags2 & flags2On) != flags2On) || (flags & flagsOff) || (flags2 & flags2Off)) continue; // Keep "bright" objects only if (nSigma > 0) if (itr->counts_model(AR_R).err() < 0 || itr->counts_model(AR_R).value() < nSigma * itr->counts_model(AR_R).err()) continue; // Calibrate object and add to appropriate chain TA_CALIB_OBJ *obj = (TA_CALIB_OBJ*) shMalloc(sizeof(TA_CALIB_OBJ)); if (shChainElementAddByPos(objChain, obj, "TA_CALIB_OBJ", TAIL, AFTER) != SH_SUCCESS) { shErrStackPush("tsFieldObjectsGet: couldn't append object chain"); shFree(obj); return SH_GENERIC_ERROR; } int raw = 1; if (tsObjCalibrate(itr, obj, NULL, (const TRANS **)transPtr1, (const TRANS **)photoTransPtr1, node, incl, 'r', 0, raw, 0) != SH_SUCCESS) { shErrStackPush("tsFieldObjectsGet: couldn't calibrate object"); return SH_GENERIC_ERROR; } obj->field = fieldInfo1; obj->run = runId; obj->camCol = camCol; obj->rerun = rerun; obj->fieldId = field->field(); } // Return return SH_SUCCESS; } /*============================================================================ ** ** ** ROUTINE: tsObjectStatusSet ** ** ** Set the SET, GOOD, and OK_RUN flags as appropriate for each ** object in a field object list. Two arrays are returned containing ** those GOOD objects (STARS, GALAXIES, and UNK only) which can ** potentially match objects in an adjacent ** field (i.e., they lie in the overlap region or a buffer around the ** overlap region of width equal to the matching radius). One array is ** an ooTVArray of object handles (ooTVArray(ooHandle(arObject))), and ** the other is an array of structures containing each object's row and ** column position along with the index of the object in the first array. ** This second array is sorted by increasing row to facilitate matching ** by the calling routine. Returns the number of potentially matchable ** objects (i.e., the size of the output arrays), or -1 on error. ** ** ** **============================================================================ */ int tsObjectStatusSet(const ooHandle(arObjectList)& field /* IN: field obj list */, int rowsPerField /* IN: non-overlap rows per field */, int overlapRows /* IN: overlap rows per field */, double matchRadius /* IN: match radius (pixels) */, double nSigma /* IN: Compile matching statistics for * nSigma detections (PSF, r') only */, ooTVArray(ooHandle(arObject)) **handles /* OUT: pointer to array of object handles */, TSOBJ **array /* OUT: pointer to array of objects, * sorted in order of increasing row */, TSMATCHSTATS *stats /* MODIFIED: matching statistics */) { // Initialize the output arrays of matchable objects. int arrayIncr = 500; int arraySize = arrayIncr; if ((*handles = new ooTVArray(ooHandle(arObject))(arraySize)) == 0) { shErrStackPush("tsObjectStatusSet: couldn't allocate memory"); return -1; } if ((*array = (TSOBJ*) shMalloc(arraySize * sizeof(TSOBJ))) == NULL) { shErrStackPush("tsObjectStatusSet: couldn't allocate memory"); return -1; } int nMatchable = 0; // Offset to convert row from field coordinate to scan coordinate int offset = field->field() * rowsPerField; // Define ranges of rows in which objects could not possibly match objects // in adjacent fields. Any objects outside this region will be returned // in the arrays. double unmatchableMin = overlapRows + matchRadius; double unmatchableMax = rowsPerField - matchRadius; // We divide the overlap regions into two, and assign // half of each overlap region to the adjacent fields, and half to the // primary region of this field. We then consider an object to belong to a // given field (OK_RUN) if it falls between those limits. double okMin = overlapRows / 2.; double okMax = rowsPerField + overlapRows / 2.; // Get iterator to objects for this field ooItr(arObject) itr; if (field->objects(itr, oocUpdate) != oocSuccess) { // Error. delete *handles; shFree(*array); shErrStackPush("tsObjectStatusSet: error getting objects for field %d", field); return -1; } // For each object ... while (itr.next()) { // Set the SET flag for all objects, indicating its been through this // step. int status = AR_OBJECT_STATUS_SET; // Determine if object is GOOD. int good = taGood(itr->objc_flags(),itr->nChild()); // If object is not GOOD, then just set its flags and skip to the // next object --- processing for BAD objects stops here. if (! good) { itr->setStatus(0, status); itr->setStatus(1, status); continue; } status |= AR_OBJECT_STATUS_GOOD; // Determine whether object lies in the primary or overlap regions for // this field and set all flags. double row = itr->objc_rowc().value(); if (row >= okMin && row < okMax) status |= AR_OBJECT_STATUS_OK_RUN; itr->setStatus(0, status); itr->setStatus(1, status); // Count only stars, galaxies, and unknown arObjectType objType = itr->objc_type(); TSSTATS *matchStats; arFloatErr counts; if (objType == AR_OBJECT_TYPE_STAR) { matchStats = &stats->stars; matchStats->nObjects++; counts = itr->psfCounts(AR_R); } else if (objType == AR_OBJECT_TYPE_GALAXY) { matchStats = &stats->galaxies; matchStats->nObjects++; counts = itr->petroCounts(AR_R); } else if (objType == AR_OBJECT_TYPE_UNK) { stats->nUnknown++; counts = itr->psfCounts(AR_R); } else continue; // Is this a "bright" object? int bright = (counts.value() > nSigma * counts.err() && counts.err() > 0.); if (bright) { if (objType == AR_OBJECT_TYPE_UNK) stats->nUnknownBright++; else matchStats->nObjectsBright++; } // If the object can't possibly match another object in an adjacent field // (i.e., its outside the overlap regions, including a buffer zone // 'matchRadius' wide), then skip to the next object. if (row > unmatchableMin && row < unmatchableMax) continue; // This object can potentially match an object in an adjacent field. // Extend the array if necessary. if (nMatchable == arraySize) { arraySize += arrayIncr; if ((*handles)->resize(arraySize) != oocSuccess) { delete *handles; shFree(*array); shErrStackPush("tsObjectStatusSet: couldn't reallocate memory"); return -1; } if ((*array = (TSOBJ*) shRealloc(*array, arraySize * sizeof(TSOBJ))) == NULL) { delete *handles; shErrStackPush("tsObjectStatusSet: couldn't reallocate memory"); return -1; } } // Add object to the array of potentially matching objects. (*handles)->set(nMatchable, itr); (*array)[nMatchable].index = nMatchable; (*array)[nMatchable].bright = bright; (*array)[nMatchable].row = row + offset; (*array)[nMatchable].col = itr->objc_colc().value(); nMatchable++; } // Resize arrays to final size if (nMatchable < arraySize) { if ((*handles)->resize(nMatchable) != oocSuccess) { delete *handles; shFree(*array); shErrStackPush("tsObjectStatusSet: couldn't resize varray"); return -1; } if ((*array = (TSOBJ*) shRealloc(*array, nMatchable * sizeof(TSOBJ))) == NULL) { if (nMatchable > 0) { delete *handles; shErrStackPush("tsObjectStatusSet: couldn't reallocate memory"); return -1; } } } // Sort it by increasing row qsort((void *)(*array), nMatchable, sizeof(TSOBJ), cmpRow); // Return number of potentially matchable objects return nMatchable; } /*============================================================================ ** ** ** ROUTINE: tsObjectStatusReset ** ** ** Reset the SET, GOOD, and OK_RUN bits in the STATUS bit mask for ** each object in a field object list. 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. **

** The number of new GOOD objects, new 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 is updated. ** ** ** **============================================================================ */ RET_CODE tsObjectStatusReset(const ooHandle(arObjectList)& field /* IN: object list */, int rowsPerField /* IN: non-overlap rows per field */, int overlapRows /* IN: overlap rows per field */, double matchRadius /* IN: match radius (pixels) */, int skyVersion /* IN: sky version to reset (0=drilling, * 1=current */, int *nNewGood /* MOD: number of new GOOD objects */, int *nNewOKRun /* MOD: number of new OK_RUN objects */, int *nOldOKRun /* MOD: number of old OK_RUN objects */, int *nProblems /* MOD: number of new OK_RUN objects which * could have matched OK_RUN objects * in adjacent fields but which won't * be border resolved */) { // We divide the overlap regions into two, and assign // half of each overlap region to the adjacent fields, and half to the // primary region of this field. We then consider an object to belong to a // given field (OK_RUN) if it falls between those limits. double okMin = overlapRows / 2.; double okMax = rowsPerField + overlapRows / 2.; // Get iterator to objects for this field ooItr(arObject) itr; if (field->objects(itr, oocUpdate) != oocSuccess) { // Error. shErrStackPush("tsObjectStatusReset: error getting objects for field %d: MUST ROLLBACK", field); return SH_GENERIC_ERROR; } // For each object ... while (itr.next()) { // Determine if object is GOOD. int good = taGood(itr->objc_flags(), itr->nChild()); // We're only concerned with objects which are GOOD now but were BAD // originally. if (! good) { if (itr->status(skyVersion) & AR_OBJECT_STATUS_GOOD) { shErrStackPush("tsObjectStatusReset: old GOOD object which is now BAD: MUST ROLLBACK"); return SH_GENERIC_ERROR; } continue; } if (itr->status(skyVersion) & AR_OBJECT_STATUS_GOOD) { if (itr->status(skyVersion) & AR_OBJECT_STATUS_OK_RUN) (*nOldOKRun)++; continue; } // This object was originally BAD but is now GOOD. We need to set its // GOOD bit. If its in the primary field region, we also need to set // its OK_RUN bit. Since we're not rematching to adjacent fields, // we can't check whether a new OK_RUN object in this field matches a new // or old OK_RUN object in an adjacent field, but we can at least // count the number of such potentially problem objects. int status = AR_OBJECT_STATUS_SET | AR_OBJECT_STATUS_GOOD; (*nNewGood)++; // Determine whether object lies in the primary or overlap regions for // this field and set all flags. double row = itr->objc_rowc().value(); if (row >= okMin && row < okMax) { status |= AR_OBJECT_STATUS_OK_RUN; (*nNewOKRun)++; } itr->setStatus(skyVersion, status); // If the object could potentially be a duplicate OK_RUN object with // another OK_RUN object in an adjacent field update the could of // potential problem cases. if ((row >= okMin - matchRadius && row <= okMin + matchRadius) || (row >= okMax - matchRadius && row <= okMax + matchRadius)) (*nProblems)++; } return SH_SUCCESS; } /********************************************************************/ /* qsort helper function. Final list is ordered from smallest to largest */ static int cmpRow (const void *d1, const void *d2) { TSOBJ *f1 = (TSOBJ *)d1; TSOBJ *f2 = (TSOBJ *)d2; if (f1->row < f2->row) return -1; if (f1->row > f2->row) return 1; return 0; } /*============================================================================ ** ** ** ROUTINE: tsSelfMerge ** ** ** Merge objects in overlap region of adjacent fields in a Frames ** Pipeline Run. Objects match if their separation ** is less than or equal to the match radius. This resolves any matched ** pairs such that, for each pair, one object is set OK_RUN and the other ** aint. All objects in a pair have their DUPLICATE flags set. Matching ** statistics are compiled and returned. The input arrays ("array" and ** "previousArray") must be sorted by increasing row. ** ** ** **============================================================================ */ void tsSelfMerge(ooTVArray(ooHandle(arObject)) *handles /* MODIFIED: array of object handles for the field */, TSOBJ *array /* IN: Array of objects for the field */, ooTVArray(ooHandle(arObject)) *previousHandles /* MODIFIED: array of object handles for the previous field */, TSOBJ *previousArray /* IN: Array of objects for the previous field */, double matchRadius /* IN: Match radius (pixels) */, double rowMin /* IN: Minimum row to match */, double rowMax /* IN: Maximum row to match */, TSMATCHSTATS *stats /* MODIFIED: matching statistics */) { // Return if either field is empty, and therefore no possible matches. int size = handles->size(); int previousSize = previousHandles->size(); if (size == 0 || previousSize == 0) return; // Set the lower end of the matching range double negMatchRadius = -matchRadius; double matchRadiusSq = matchRadius * matchRadius; // Find index to start searches in previous field (just searching in upper // overlap region). Limit array to just that portion. Since current field // is being searched in the lower overlap region, we start at the beginning // of that array. int previousStart = tsObjArrayBinomialSearch(rowMin, previousArray, previousSize) + 1; if (previousStart == previousSize) return; // no objects in overlap region // For each object in this field ... for (int i = 0; i < size; i++) { // Both arrays are sorted by increasing row, so we search on row. double row = array[i].row; // If beyond max row then we have left the lower overlap area and are // done. if (row > rowMax) break; // Increment low index from previous search until we hit mininum row // for this search. while (previousStart < previousSize) { if (previousArray[previousStart].row >= row - matchRadius) break; previousStart++; } // Skip if all objects in the previous array are at a smaller row. // In fact, since this array is ordered by increasing row, we can // break out of the loop since there can be no more matches. if (previousStart == previousSize) break; // Now check each object until we're beyond the match radius or we // run out of the array. int closestIndex = -1; double closestDistance = 1000. * matchRadiusSq; double rowHi = row + matchRadius; double col = array[i].col; for (int index = previousStart; index < previousSize; index++) { // Break out if we're beyond the match radius if (previousArray[index].row > rowHi) break; double colDiff = col - previousArray[index].col; if (colDiff > matchRadius || colDiff < negMatchRadius) continue; double rowDiff = row - previousArray[index].row; double distance = rowDiff * rowDiff + colDiff * colDiff; if (distance <= matchRadiusSq) { // This matches. Save it if its the closest matching one // thus far. if (distance < closestDistance) { closestIndex = index; closestDistance = distance; } } } // Do we have a match? if (closestIndex != -1) { // Yes. Link the two objects in the database. (*handles)[array[i].index]-> addDuplicates((*previousHandles)[previousArray[closestIndex].index]); // Resolve the pair, so that only one object in the pair is set to // OK_RUN. Since an object in this field can have at most one match, // while objects in the previous field can have more than one match, // we keep the OK_RUNedness of the object in the previous field and // adjust the object in this field accordingly. int status = (*handles)[array[i].index]->status(0); int previousStatus = (*previousHandles)[previousArray[closestIndex].index]->status(0); if ((status & AR_OBJECT_STATUS_OK_RUN) == (previousStatus & AR_OBJECT_STATUS_OK_RUN)) { if (previousStatus & AR_OBJECT_STATUS_OK_RUN) status &= ~AR_OBJECT_STATUS_OK_RUN; else status |= AR_OBJECT_STATUS_OK_RUN; } (*handles)[array[i].index]-> setStatus(0, status | AR_OBJECT_STATUS_DUPLICATE); (*handles)[array[i].index]-> setStatus(1, status | AR_OBJECT_STATUS_DUPLICATE); (*previousHandles)[previousArray[closestIndex].index]-> setStatus(0, previousStatus | AR_OBJECT_STATUS_DUPLICATE); (*previousHandles)[previousArray[closestIndex].index]-> setStatus(1, previousStatus | AR_OBJECT_STATUS_DUPLICATE); // Update matching statistics arObjectType objType = (*handles)[array[i].index]->objc_type(); TSSTATS *matchStats; if (objType == AR_OBJECT_TYPE_STAR) matchStats = &stats->stars; else if (objType == AR_OBJECT_TYPE_GALAXY) matchStats = &stats->galaxies; else continue; matchStats->nMatches++; if (objType != (*previousHandles)[previousArray[closestIndex].index] ->objc_type()) { stats->nMismatches++; if (array[i].bright) stats->nMismatchesBright++; } // Add to statistics only if object is "bright" (at least an // nSigma detection). if (array[i].bright) { matchStats->nMatchesBright++; double rowDiff = row - previousArray[closestIndex].row; double colDiff = col - previousArray[closestIndex].col; matchStats->row.av += rowDiff; matchStats->row.rms += rowDiff * rowDiff; matchStats->row.nPairs++; matchStats->col.av += colDiff; matchStats->col.rms += colDiff * colDiff; matchStats->col.nPairs++; arFloatErr model1, model2, bModel1, bModel2; for (int k = 0; k < nFilters; k++) { arFloatErr counts1, counts2; if (objType == AR_OBJECT_TYPE_STAR) { // Use PSF magnitudes for stars counts1 = (*handles)[array[i].index]->psfCounts((arFilter)k); counts2 = (*previousHandles)[previousArray[closestIndex].index] ->psfCounts((arFilter)k); } else { // Use petrosian magnitudes for galaxies counts1 = (*handles)[array[i].index]->petroCounts((arFilter)k); counts2 = (*previousHandles)[previousArray[closestIndex].index] ->petroCounts((arFilter)k); } if (counts1.err() >= 0 && counts2.err() >= 0) { double magDiff = luptitude(counts1.value()) - luptitude(counts2.value()); matchStats->mag[k].av += magDiff; matchStats->mag[k].rms += magDiff * magDiff; matchStats->mag[k].nPairs++; } // Update color stats model1 = (*handles)[array[i].index]->counts_model((arFilter)k); model2 = (*previousHandles)[previousArray[closestIndex].index] ->counts_model((arFilter)k); if (k > 0) if (model1.err() > 0 && model2.err() > 0 && bModel1.err() > 0 && bModel2.err() > 0) { double magDiff = (luptitude(bModel1.value()) - luptitude(model1.value())) - (luptitude(bModel2.value()) - luptitude(model2.value())); matchStats->color[k-1].av += magDiff; matchStats->color[k-1].rms += magDiff * magDiff; matchStats->color[k-1].nPairs++; } bModel1 = model1; bModel2 = model2; } } } } // Final star statistics if (stats->stars.row.nPairs > 0) { stats->stars.row.av /= stats->stars.row.nPairs; double meanSq = stats->stars.row.rms / stats->stars.row.nPairs - stats->stars.row.av * stats->stars.row.av; stats->stars.row.rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } if (stats->stars.col.nPairs > 0) { stats->stars.col.av /= stats->stars.col.nPairs; double meanSq = stats->stars.col.rms / stats->stars.col.nPairs - stats->stars.col.av * stats->stars.col.av; stats->stars.col.rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } for (i = 0; i < nFilters; i++) if (stats->stars.mag[i].nPairs > 0) { stats->stars.mag[i].av /= stats->stars.mag[i].nPairs; double meanSq = stats->stars.mag[i].rms / stats->stars.mag[i].nPairs - stats->stars.mag[i].av * stats->stars.mag[i].av; stats->stars.mag[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } for (i = 0; i < nColors; i++) if (stats->stars.color[i].nPairs > 0) { stats->stars.color[i].av /= stats->stars.color[i].nPairs; double meanSq = stats->stars.color[i].rms / stats->stars.color[i].nPairs - stats->stars.color[i].av * stats->stars.color[i].av; stats->stars.color[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } // Final galaxy statistics if (stats->galaxies.row.nPairs > 0) { stats->galaxies.row.av /= stats->galaxies.row.nPairs; double meanSq = stats->galaxies.row.rms / stats->galaxies.row.nPairs - stats->galaxies.row.av * stats->galaxies.row.av; stats->galaxies.row.rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } if (stats->galaxies.col.nPairs > 0) { stats->galaxies.col.av /= stats->galaxies.col.nPairs; double meanSq = stats->galaxies.col.rms / stats->galaxies.col.nPairs - stats->galaxies.col.av * stats->galaxies.col.av; stats->galaxies.col.rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } for (i = 0; i < nFilters; i++) if (stats->galaxies.mag[i].nPairs > 0) { stats->galaxies.mag[i].av /= stats->galaxies.mag[i].nPairs; double meanSq = stats->galaxies.mag[i].rms / stats->galaxies.mag[i].nPairs - stats->galaxies.mag[i].av * stats->galaxies.mag[i].av; stats->galaxies.mag[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } for (i = 0; i < nColors; i++) if (stats->galaxies.color[i].nPairs > 0) { stats->galaxies.color[i].av /= stats->galaxies.color[i].nPairs; double meanSq = stats->galaxies.color[i].rms / stats->galaxies.color[i].nPairs - stats->galaxies.color[i].av * stats->galaxies.color[i].av; stats->galaxies.color[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.); } // Return. return; } /*============================================================================ ** ** ** ROUTINE: tsRunToCalibArray ** ** ** Create an array of calibrated objects (TSCALIBOBJ) for all OK_RUN ** objects (STARS, GALAXIES, and UNK only) in a Frames Pipeline Run. The ** array will be sorted in increasing mu order. ** DCR corrections are applied using cosmic colors, since ** a final photometric calibration can not be guaranteed to exist for ** the run at this time. ** Returns the number of objects read, or -1 on error. ** ** ** **============================================================================ */ int tsRunToCalibArray(const ooHandle(arFramesPipeRun)& run /* IN: Frames Pipeline run */, const ooHandle(arAstromCalib)& astrom /* IN: astrometric calibration */, ooTVArray(ooHandle(arObject)) **handles /* OUT: pointer to array of object handles */, TSCALIBOBJ **array /* OUT: pointer to array of calibrated objects, sorted in * order of increasing mu */, double *minNu /* OUT: minimum nu in array (radians) */, double *maxNu /* OUT: maximum nu in array (radians) */) { // Initialize the array of calibrated objects. int arrayIncr = 50000; int arraySize = arrayIncr; if ((*handles = new ooTVArray(ooHandle(arObject))(arraySize)) == 0) { shErrStackPush("tsRunToCalibArray: couldn't allocate memory"); return NULL; } if ((*array = (TSCALIBOBJ*) shMalloc(arraySize * sizeof(TSCALIBOBJ))) == NULL) { shErrStackPush("tsRunToCalibArray: couldn't allocate memory"); return NULL; } int nRead = 0; // Fetch the range of fields uint16 field0 = run->field0(); uint16 nFields = run->nFields(); uint8 camCol = run->camCol(); int first = 1; // 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 ccd = 10 * camRow + camCol; // For each field ... for (int32 field = field0; field < field0 + nFields; field++) { printf("%4d", field); fflush(stdout); // Get the astrometric transformation for this field TRANS trans = astrom->astrotoolsTrans(ccd, field); // Get iterator to objects for this field ooHandle(arObjectList) objectList = run->objectList(field); ooItr(arObject) itr; if (objectList->objects(itr, oocUpdate) != oocSuccess) { // Error. delete *handles; shFree(*array); shErrStackPush("tsRunToCalibArray: error getting objects for field %d", field); return -1; } // For each object ... while (itr.next()) { // We will match only OK_RUN objects. if (! (itr->status(0) & AR_OBJECT_STATUS_OK_RUN)) continue; // Only match stars, galaxies, and unknowns arObjectType objType = itr->objc_type(); if (objType != AR_OBJECT_TYPE_STAR && objType != AR_OBJECT_TYPE_GALAXY && objType != AR_OBJECT_TYPE_UNK) continue; // Extend the array if necessary if (nRead == arraySize) { arraySize += arrayIncr; if ((*handles)->resize(arraySize) != oocSuccess) { delete *handles; shFree(*array); shErrStackPush("tsRunToCalibArray: couldn't reallocate memory"); return -1; } if ((*array = (TSCALIBOBJ*) shRealloc(*array, arraySize * sizeof(TSCALIBOBJ))) == NULL) { delete *handles; shErrStackPush("tsRunToCalibArray: couldn't reallocate memory"); return -1; } } // Append calibrated object to the array (*handles)->set(nRead, itr); arSphericalCoord coord = itr->gcCoord(trans); (*array)[nRead].index = nRead; (*array)[nRead].mu = coord.longitude(); (*array)[nRead].nu = coord.latitude(); // Update minimum and maximum nu in the array if (first) { *minNu = (*array)[nRead].nu; *maxNu = (*array)[nRead].nu; first = 0; } else { if ((*array)[nRead].nu < *minNu) *minNu = (*array)[nRead].nu; if ((*array)[nRead].nu > *maxNu) *maxNu = (*array)[nRead].nu; } // Increment the total number of objects nRead++; } } printf("\n"); // Resize the array if (nRead < arraySize) { if ((*handles)->resize(nRead) != oocSuccess) { delete *handles; shFree(*array); shErrStackPush("tsRunToCalibArray: couldn't reallocate memory"); return -1; } if ((*array = (TSCALIBOBJ*) shRealloc(*array, nRead * sizeof(TSCALIBOBJ))) == NULL) { delete *handles; shErrStackPush("tsRunToCalibArray: couldn't reallocate memory"); return -1; } } // Sort it by increasing mu qsort((void *)(*array), nRead, sizeof(TSCALIBOBJ), cmpMu); // Return number of objects read return nRead; } /********************************************************************/ /* qsort helper function. Final list is ordered from smallest to largest*/ static int cmpMu (const void *d1, const void *d2) { TSCALIBOBJ *f1 = (TSCALIBOBJ *)d1; TSCALIBOBJ *f2 = (TSCALIBOBJ *)d2; if (f1->mu < f2->mu) return -1; if (f1->mu > f2->mu) return 1; return 0; } /*============================================================================ ** ** ** ROUTINE: tsRunMergeOne ** ** ** Merge one Frames Pipeline Run against another. Return SH_SUCCESS ** if successful, else an error code. Only OK_RUN objects of type ** STAR, GALAXY, or UNK are matched. **

** Assumes that nu is sufficiently small that cos(nu) ~ 1 and can thus ** be ignored. DCR corrections are applied using cosmic colors, since ** a final photometric calibration can not be guaranteed to exist for ** the run at this time. ** ** ** **============================================================================ */ RET_CODE tsRunMergeOne(ooTVArray(ooHandle(arObject)) *handles /* MODIFIED: array of object handles */, TSCALIBOBJ *array /* MODIFIED: Array of calibrated objects */, int nObjects /* IN: Number of objects to merge */, double minNu /* IN: minimum nu in array (degrees) */, double maxNu /* IN: maximum nu in array (degrees) */, double node /* IN: ascending node of GC (degrees) */, double incl /* IN: inclination of GC (degrees) */, double equinox /* IN: equinox of GC coords (years) */, ooHandle(arFramesPipeRun)& other /* IN: Frames Pipe Run to merge against */, double matchRadius /* IN: Match radius (arcsecs) */) { // Fetch the range of fields and calibrations for the other run uint16 field0 = other->field0(); uint16 nFields = other->nFields(); uint8 camCol = other->camCol(); ooHandle(arAstromCalib) astromH = other->mergeAstrometricCalibration(); if (fabs(astromH->equinox() - equinox) > 0.000001) { shErrStackPush("tsRunMergeOne: mismatched equinoxes"); return SH_GENERIC_ERROR; } double node2 = astromH->node(); double incl2 = astromH->incl(); // 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 ccd = 10 * camRow + camCol; // Get number of non-overlap rows and columns on a frame. int nRows = other->imagingRun()->program()->ccd(camRow,other->camCol())->nRows(); // HARDWIRED: number of columns per frame. int nCols = 2048; // HARDWIRED: number of overlap rows per frame. int overlapRows = 128; // Want full size of the frame nRows += overlapRows; // Convert match radius to degrees matchRadius *= at_asec2Deg; double negMatchRadius = -matchRadius; double matchRadiusSq = matchRadius * matchRadius; // Range in great circle longitude covered by the array (degrees) double minMu = array[0].mu; double maxMu = array[nObjects-1].mu; // Some parked scans have oddly defined great circles in which the scans // can go through 0/360. Thus, to match-up, we'll need to force the "other" // scan into the same range of mu as "this" scan. Thus, lets determine a // good range of mu for "this" scan (just subtract 10 degrees from the // minimum mu for the low mu in the range). double muRangeLo = minMu - 10.; double muRangeHi = muRangeLo + 360.; // For each field ... printf("run %d, camCol %d, rerun %d, fields %d to %d\n", other->run(), other->camCol(), other->rerun(), field0, field0+nFields-1); fflush(stdout); for (int32 field = field0; field < field0 + nFields; field++) { // Get astrometric transformation for this field TRANS trans = astromH->astrotoolsTrans(ccd, field); // Get great circle coordinates of the field center (degrees). double centerMu = 0, centerNu = 0; atTransApply(&trans, 'r', nRows / 2., 0., nCols / 2., 0., NULL, NULL, ¢erMu, NULL, ¢erNu,NULL); // If first field, measure the field radius (degrees). double fieldRadius; if (field == field0) { double cornerMu = 0, cornerNu = 0; atTransApply(&trans, 'r', 0., 0., 0., 0., NULL, NULL, &cornerMu, NULL, &cornerNu, NULL); fieldRadius = arSphericalCoord(centerMu, centerNu). distance(arSphericalCoord(cornerMu, cornerNu)) * at_asec2Deg; } // Convert center coordinates to target run great circle coordinates // (degrees). double ra, dec; atGCToEq(centerMu, centerNu, &ra, &dec, node2, incl2); atEqToGC(ra, dec, ¢erMu, ¢erNu, node, incl); atBound(¢erMu, muRangeLo, muRangeHi); // Skip this field if it lies outside range of mu or nu spanned by the // target run. if (centerMu + fieldRadius < minMu || centerMu - fieldRadius > maxMu) continue; if (centerNu + fieldRadius < minNu || centerNu - fieldRadius > maxNu) continue; printf("%4d", field); fflush(stdout); // Get iterator to objects for this field ooHandle(arObjectList) objectList = other->objectList(field); ooItr(arObject) itr; if (objectList->objects(itr, oocUpdate) != oocSuccess) { shErrStackPush("tsRunMergeOne: error getting objects for field %d", field); return SH_GENERIC_ERROR; } // For each object ... int first = 1; while (itr.next()) { // Only match OK_RUN objects. if (! (itr->status(0) & AR_OBJECT_STATUS_OK_RUN)) continue; // Only match stars, galaxies, and unknowns arObjectType objType = itr->objc_type(); if (objType != AR_OBJECT_TYPE_STAR && objType != AR_OBJECT_TYPE_GALAXY && objType != AR_OBJECT_TYPE_UNK) continue; // Fetch calibrated great circle coordinates (degrees) arSphericalCoord coord = itr->gcCoord(trans); double mu = coord.longitude(); double nu = coord.latitude(); // Convert to target run great circle coordinates (degrees) atGCToEq(mu, nu, &ra, &dec, node2, incl2); atEqToGC(ra, dec, &mu, &nu, node, incl); atBound(&mu, muRangeLo, muRangeHi); // Look for match in calibrated array. We'll search based on the // object position minus the match radius. int index; if (first) { // Do a binomial search for the first object. The binomial // search returns an array index which is set to one below the // lowest matching or greater than array element. index = tsCalibObjArrayBinomialSearch(mu-matchRadius, array, nObjects) + 1; first = 0; } else { // Search with correlated values around the last index. index = tsCalibObjArrayCorrelatedSearch(mu-matchRadius, array, nObjects, index-1) + 1; } // Skip if all objects in the array are at a smaller mu. if (index == nObjects) continue; // Now check each object until we're beyond the match radius or we // run out of the array. int closestIndex = -1; double closestDistance = 1000. * matchRadiusSq; double muHi = mu + matchRadius; for (; index < nObjects; index++) { // Break out if we're beyond the match radius if (array[index].mu > muHi) break; double nuDiff = nu - array[index].nu; if (nuDiff > matchRadius || nuDiff < negMatchRadius) continue; double muDiff = mu - array[index].mu; double distance = muDiff * muDiff + nuDiff * nuDiff; if (distance <= matchRadiusSq) { // This matches. Save it if its the closest matching one // thus far. if (distance < closestDistance) { closestIndex = index; closestDistance = distance; } } } // If we had a match, link the two objects. if (closestIndex != -1) (*handles)[array[closestIndex].index]->addMatches(itr); } } printf("\n"); fflush(stdout); // Return. return SH_SUCCESS; } // Binomial search by mu on an array of TSCALIBOBJ. Returns the index such // that the input mu is in the range mu(index) to mu(index+1); // the input mu is definitely greater than mu(index), but could be // equal to mu(index+1). int tsCalibObjArrayBinomialSearch(double mu, TSCALIBOBJ *array, int nObjects) { int lo = -1; int hi = nObjects; while (hi - lo > 1) { int mid = (hi + lo) / 2; if (mu > array[mid].mu) lo = mid; else hi = mid; } return lo; } // Binomial search by row on an array of TSOBJ. Returns the index such // that the input row is in the range row(index) to row(index+1); // the input row is definitely greater than row(index), but could be // equal to row(index+1). int tsObjArrayBinomialSearch(double row, TSOBJ *array, int nObjects) { int lo = -1; int hi = nObjects; while (hi - lo > 1) { int mid = (hi + lo) / 2; if (row > array[mid].row) lo = mid; else hi = mid; } return lo; } // Correlated search by mu on an array of TSCALIBOBJ. Returns the index such // that the input mu is in the range mu(index) to mu(index+1); // the input mu is definitely greater than mu(index), but could be // equal to mu(index+1). int tsCalibObjArrayCorrelatedSearch(double mu, TSCALIBOBJ *array, int nObjects, int guess) { if (guess <= -1) { if (mu <= array[0].mu) return -1; else guess = 0; } else if (guess >= nObjects) { if (mu > array[nObjects - 1].mu) return nObjects - 1; else guess = nObjects - 1; } int lo, hi; int incr = 1; if (mu > array[guess].mu) { if (guess == nObjects - 1) return guess; lo = guess; hi = guess + incr; while (mu > array[hi].mu) { lo = hi; incr += incr; hi += incr; if (hi >= nObjects) { hi = nObjects; break; } } } else { if (guess == 0) return -1; hi = guess; lo = guess - incr; while (mu <= array[lo].mu) { hi = lo; incr += incr; lo -= incr; if (lo < 0) { lo = -1; break; } } } while (hi - lo > 1) { int mid = (hi + lo) / 2; if (mu > array[mid].mu) lo = mid; else hi = mid; } return lo; }