Index: /branches/eam_branches/ipp-20140717/Ohana/src/relphot/include/relphot.h =================================================================== --- /branches/eam_branches/ipp-20140717/Ohana/src/relphot/include/relphot.h (revision 37112) +++ /branches/eam_branches/ipp-20140717/Ohana/src/relphot/include/relphot.h (revision 37113) @@ -440,4 +440,5 @@ int MatchImageName (off_t meas, int cat, char *name); int MatchImageSkycellID (off_t meas, int cat, int myTessID, int myProjectionID, int mySkycellID); +int FindImageSkycellID (off_t meas, int cat, int *myTessID, int *myProjectionID, int *mySkycellID); int load_tess (char *treefile); Index: /branches/eam_branches/ipp-20140717/Ohana/src/relphot/src/ImageOps.c =================================================================== --- /branches/eam_branches/ipp-20140717/Ohana/src/relphot/src/ImageOps.c (revision 37112) +++ /branches/eam_branches/ipp-20140717/Ohana/src/relphot/src/ImageOps.c (revision 37113) @@ -491,4 +491,27 @@ if (projectID[i] != myProjectionID) return FALSE; if (skycellID[i] != mySkycellID) return FALSE; + + return TRUE; +} + +// returns image.Mcal - ff(x,y) +int FindImageSkycellID (off_t meas, int cat, int *myTessID, int *myProjectionID, int *mySkycellID) { + + off_t i; + + if (!MeasureToImage) return FALSE; + + if (meas < 0) return FALSE; + + i = MeasureToImage[cat][meas]; + if (i == -1) return FALSE; + + if (tessID[i] == -1) return FALSE; + if (projectID[i] == -1) return FALSE; + if (skycellID[i] == -1) return FALSE; + + *myTessID = tessID[i]; + *myProjectionID = projectID[i]; + *mySkycellID = skycellID[i]; return TRUE; Index: /branches/eam_branches/ipp-20140717/Ohana/src/relphot/src/StarOps.c =================================================================== --- /branches/eam_branches/ipp-20140717/Ohana/src/relphot/src/StarOps.c (revision 37112) +++ /branches/eam_branches/ipp-20140717/Ohana/src/relphot/src/StarOps.c (revision 37113) @@ -1086,4 +1086,10 @@ if (ecode != thisCode) { continue; } + // SKIP gpc1 stack data + if (isGPC1stack(catalog[i].measureT[m].photcode)) continue; + + // SKIP gpc1 forced-warp data + if (isGPC1warp(catalog[i].measureT[m].photcode)) continue; + // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers Mcal = getMcal (m, i, flatcorr, catalog); @@ -1131,4 +1137,10 @@ if (ecode != thisCode) { continue; } + // SKIP gpc1 stack data + if (isGPC1stack(catalog[i].measureT[m].photcode)) continue; + + // SKIP gpc1 forced-warp data + if (isGPC1warp(catalog[i].measureT[m].photcode)) continue; + // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers Mcal = getMcal (m, i, flatcorr, catalog); Index: /branches/eam_branches/ipp-20140717/Ohana/src/relphot/src/bcatalog.c =================================================================== --- /branches/eam_branches/ipp-20140717/Ohana/src/relphot/src/bcatalog.c (revision 37112) +++ /branches/eam_branches/ipp-20140717/Ohana/src/relphot/src/bcatalog.c (revision 37113) @@ -78,5 +78,5 @@ int Nsec = GetPhotcodeNsec(thisCode); - fprintf (stderr, "really use dvo_secfilt_init?"); + // fprintf (stderr, "really use dvo_secfilt_init?"); dvo_secfilt_init (&subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec]); Index: /branches/eam_branches/ipp-20140717/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- /branches/eam_branches/ipp-20140717/Ohana/src/relphot/src/setMrelCatalog.c (revision 37112) +++ /branches/eam_branches/ipp-20140717/Ohana/src/relphot/src/setMrelCatalog.c (revision 37113) @@ -441,9 +441,13 @@ // only apply Stack operation on setMrelFinal in first pass +// this function has 3 goals, not to be confused: +// 1) find and mark the PRIMARY detections (regardless of the quality of the detection) +// 2) select the BEST detections per filter (regardless of PRIMARY) +// 3) apply the zero point and AB->Jy transformations int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure) { off_t k, ID; - float Msys = 0, Mcal= 0, Mmos = 0, Mgrid = 0; + float Mcal= 0, Mmos = 0, Mgrid = 0, Finst = 0; // set the primary projection cell and skycell for this coordinate @@ -467,12 +471,8 @@ int haveStack = FALSE; - // need to find the measurement closest to the center of its skycell, as well as the - // closest for the subset of primary projection cells - - float stackCenterOffsetMin = 1e9; - off_t stackCenterMeasureMin = -1; - - float stackPrimaryOffsetMin = 1e9; - off_t stackPrimaryMeasureMin = -1; + float psfQFbest = 0.0; + + off_t stackBestMeasure = -1; + off_t stackPrimaryMeasure = -1; int isBad = FALSE; @@ -481,4 +481,7 @@ off_t meas = measureOffset; for (k = 0; k < average[0].Nmeasure; k++, meas++) { + + // only examine gpc1 stack data + if (!isGPC1stack(measure[k].photcode)) continue; // skip measurements that do not match the current photcode @@ -489,13 +492,10 @@ // clear this bit for all measurements measure[k].dbFlags &= ~ID_MEAS_STACK_PRIMARY; - - // only examine gpc1 stack data - if (!isGPC1stack(measure[k].photcode)) continue; + measure[k].dbFlags &= ~ID_MEAS_STACK_PHOT_SRC; haveStack = TRUE; haveStackObject = TRUE; - if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad); - + // match measurement to its image (this is just a check, right?) if (getImageEntry (meas, cat) < 0) { fprintf (stderr, "*** warning : this should be impossible\n"); @@ -503,4 +503,19 @@ continue; } + + // find the PRIMARY measurement + + // if we request the primary (USE_TREE_FOR_PRIMARY), this is true if the measurement is from the + // primary skycell for this position + if (MatchImageSkycellID (meas, cat, tessID, projectID, skycellID)) { + stackPrimaryMeasure = k; + measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY; + secfilt[Nsec].stackPrmryOff = meas; + } + + // now choose the BEST measurements (may also be PRIMARY) + + if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad); + // measurements without an image are either external reference photometry or // data for which the associated image has not been loaded (probably because of @@ -514,26 +529,11 @@ if (isnan(Mgrid)) SKIP_THIS_MEAS_STACK(Ngrid); - Msys = PhotSys (&measure[k], &average[0], &secfilt[0], MAG_CLASS_PSF); - if (isnan(Msys)) SKIP_THIS_MEAS_STACK(Nsys); - - unsigned int stackImageID; - - // which stack image should we use for the mean value? - // if we request the primary (USE_TREE_FOR_PRIMARY), then find the min distances for data from the primary cell - if (MatchImageSkycellID (meas, cat, tessID, projectID, skycellID)) { - float stackPrimaryOffset = getCenterOffset (meas, cat, &measure[k], &stackImageID); - if (stackPrimaryOffset < stackPrimaryOffsetMin) { - stackPrimaryOffsetMin = stackPrimaryOffset; - // stackPrimaryIDmin = stackImageID; - stackPrimaryMeasureMin = k; - } - } - - // get the center distance for the generic case: - float stackCenterOffset = getCenterOffset (meas, cat, &measure[k], &stackImageID); - if (stackCenterOffset < stackCenterOffsetMin) { - stackCenterOffsetMin = stackCenterOffset; - // stackCenterIDmin = stackImageID; - stackCenterMeasureMin = k; + // NOTE: negative and insignificant vlues values are allowed, but not NAN flux values + Finst = PhotFluxInst (&measure[k], MAG_CLASS_PSF); + if (isnan(Finst)) SKIP_THIS_MEAS_STACK(Ninst); + + if (measure[k].psfQFperf > psfQFbest) { + psfQFbest = measure[k].psfQFperf; + stackBestMeasure = k; } @@ -558,7 +558,6 @@ // measurements which are bad will not have a valid stack entry and are skipped - k = (stackPrimaryMeasureMin >= 0) ? stackPrimaryMeasureMin : stackCenterMeasureMin; + k = (stackBestMeasure >= 0) ? stackBestMeasure : stackPrimaryMeasure; if (k < 0) continue; - // we are now populating stackDetectID not stack Image ID in secfilt @@ -609,13 +608,14 @@ secfilt[Nsec].dFapStk = zpFactor * measure[k].dFluxAp; + // Jy to AB mags secfilt[Nsec].MpsfStk = (measure[k].FluxPSF > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FpsfStk) : NAN; secfilt[Nsec].MkronStk = (measure[k].FluxKron > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FkronStk) : NAN; secfilt[Nsec].MapStk = (measure[k].FluxAp > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FapStk) : NAN; - secfilt[Nsec].stackDetectID = ID; - + secfilt[Nsec].stackBestOff = k + measureOffset; + + // this is the measurement used by secfilt[] measure[k].dbFlags |= ID_MEAS_STACK_PHOT_SRC; - if (stackPrimaryMeasureMin >= 0) { - measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY; + if (k == stackPrimaryMeasure) { secfilt[Nsec].flags |= ID_SECF_STACK_PRIMARY; } @@ -646,5 +646,5 @@ off_t k; - float Fsys = 0, Mcal= 0; + float Mcal= 0; // we are measuring means for 3 types of FLUXes: psf, ap, kron. I am using the psf mag @@ -668,5 +668,5 @@ // option for a test print - if (FALSE && (average[0].objID == 0x7146) && (average[0].catID == 0x49d8)) { + if ((average[0].objID == 0xd6e) && (average[0].catID == 0x4984)) { fprintf (stderr, "test obj\n"); print_measure_set_alt (average, secfilt, measure); @@ -684,4 +684,9 @@ if (found[Nsec]) continue; + // I need to restrict the measurements coming only from the same skycell as the secfilt values + // if we cannot make that association, skip this set of warps + int stkTessID, stkProjID, stkSkycellID; + if (!FindImageSkycellID (secfilt[Nsec].stackBestOff, cat, &stkTessID, &stkProjID, &stkSkycellID)) continue; + off_t meas = measureOffset; @@ -698,4 +703,8 @@ if (!code) continue; if (code->equiv != thisCode) { continue; } + + measure[k].dbFlags &= ~ID_MEAS_WARP_USED; + + if (!MatchImageSkycellID (meas, cat, stkTessID, stkProjID, stkSkycellID)) continue; if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); @@ -724,9 +733,18 @@ // are small (and not measured) float Fpsf = PhotFluxCat (&measure[k], MAG_CLASS_PSF); - if (isnan(Fsys)) SKIP_THIS_MEAS(Nsys); + if (isnan(Fpsf)) SKIP_THIS_MEAS(Nsys); // if (Msys < 0.0) SKIP_THIS_MEAS(Nsys); // if (Msys > 30.0) SKIP_THIS_MEAS(Nsys); float dFpsf = PhotFluxCatErr (&measure[k], MAG_CLASS_PSF); + if (isnan(dFpsf)) SKIP_THIS_MEAS (Nsys); + + int isBad = (measure[k].psfQF < 0.85); + isBad |= isnan(measure[k].psfQF); + if (isBad) SKIP_THIS_MEAS (Nsys); + + int isPoor = (measure[k].psfQFperf < 0.85); + if ((pass == 0) && isPoor) SKIP_THIS_MEAS (Nsys); + dFpsf = MAX (dFpsf, MIN_ERROR*Fpsf); // MIN_ERROR is a fractional error Fpsflist[Npsf] = Fpsf * Fcal; @@ -734,4 +752,6 @@ wpsflist[Npsf] = 1.0; Npsf ++; + + measure[k].dbFlags |= ID_MEAS_WARP_USED; float Fap = PhotFluxCat (&measure[k], MAG_CLASS_APER);