Index: branches/eam_branches/ohana.20170822/src/libdvo/include/dvo.h =================================================================== --- branches/eam_branches/ohana.20170822/src/libdvo/include/dvo.h (revision 40248) +++ branches/eam_branches/ohana.20170822/src/libdvo/include/dvo.h (revision 40249) @@ -248,5 +248,5 @@ ID_SECF_CHIP_FLAGS = 0x01003f1f, // all chip-related bits (used to reset the correct bits only) - ID_SECF_STACK_FLAGS = 0x00004020, // all stack-related bits ( + ID_SECF_STACK_FLAGS = 0x0003c020, // all stack-related bits ( } DVOSecfiltFlags; Index: branches/eam_branches/ohana.20170822/src/relphot/src/setMrelCatalog.c =================================================================== --- branches/eam_branches/ohana.20170822/src/relphot/src/setMrelCatalog.c (revision 40248) +++ branches/eam_branches/ohana.20170822/src/relphot/src/setMrelCatalog.c (revision 40249) @@ -565,6 +565,4 @@ int haveStack = FALSE; - off_t stackPrimaryMeasure = -1; - int isBad = FALSE; int isSuspect = FALSE; @@ -572,4 +570,8 @@ int NstackMeas = 0; // number of stack measurements for this photcode int NstackDet = 0; // number of stack detections for this photcode (not forced) + + int primaryLevelMax = 0; + off_t primaryEntryMax = -1; + float primaryValueMax = 0.0; int bestLevelMax = 0; @@ -579,9 +581,21 @@ int Nprimary = 0; + // reset all stack-related values for this secfilt: secfilt[Nsec].stackBestOff = -1; secfilt[Nsec].stackPrmryOff = -1; - - off_t meas = measureOffset; - for (off_t k = 0; k < Nmeasure; k++, meas++) { + secfilt[Nsec].FpsfStk = NAN; + secfilt[Nsec].dFpsfStk = NAN; + secfilt[Nsec].FkronStk = NAN; + secfilt[Nsec].dFkronStk = NAN; + secfilt[Nsec].FapStk = NAN; + secfilt[Nsec].dFapStk = NAN; + secfilt[Nsec].MpsfStk = NAN; + secfilt[Nsec].MkronStk = NAN; + secfilt[Nsec].MapStk = NAN; + secfilt[Nsec].Nstack = 0; + secfilt[Nsec].NstackDet = 0; + secfilt[Nsec].flags &= ~ID_SECF_STACK_FLAGS; // reset the stack flags + + for (off_t k = 0; k < Nmeasure; k++) { // only examine gpc1 stack data @@ -594,5 +608,5 @@ NstackMeas ++; - if ((measure[k].photFlags2 & 0x00000004) == 0) NstackDet ++; + if ((measure[k].photFlags2 & 0x00000004) == 0) NstackDet ++; // 0x04 = PM_SOURCE_MODE2_MATCHED // clear this bit for all measurements @@ -603,5 +617,5 @@ haveStackObject = TRUE; - XXX ??? if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad); + // if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad); int isPrimary = FALSE; @@ -610,6 +624,7 @@ // 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)) { + // is from the primary skycell for this position. (note that MatchImageSkycellID + // requires the entry in the full Measure table, not just this object) + if (MatchImageSkycellID (measureOffset + k, cat, tessID, projectID, skycellID)) { measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY; isPrimary = TRUE; @@ -617,5 +632,7 @@ } - int psfQFperfAboveLimit = isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > 0.95); + // soften the error floor (can dM be 0.0?) + int SNvalue = isfinite(measure[k].dM) ? 1.0 / hypot (measure[k].dM, MIN_ERROR) : NAN; + int psfQFperfAboveLimit = isfinite(SNvalue) && isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > 0.95); // ** determine the BEST level @@ -627,6 +644,5 @@ // here is the rule for choosing the best value: - // XXX what if dM = 0.0 or NAN? (treat as bestLevel 1 or 2?) - bestValue = (bestLevel > 2) ? 1.0 / measure[k].dM : measure[k].psfQFperf; + float bestValue = (bestLevel > 2) ? SNvalue : measure[k].psfQFperf; // if we have not reached this level before, set the new level @@ -649,8 +665,7 @@ // here is the rule for choosing the PRIMARY value: - // XXX what if dM = 0.0 or NAN? (treat as bestLevel 1 or 2?) - primaryValue = 0.0; + float primaryValue = 0.0; if (primaryLevel == 1) { primaryValue = measure[k].psfQFperf; } - if (primaryLevel == 2) { primaryValue = 1.0 / measure[k].dM; } + if (primaryLevel == 2) { primaryValue = SNvalue; } // if we have not reached this level before, set the new level @@ -672,8 +687,5 @@ // now we have a BEST and a PRIMARY entry (may be the same entry) - secfilt[Nsec].stackPrmryOff = (primaryEntryMax == -1) ? -1 : primaryEntryMax + measureOffset; - myAssert (secfilt[Nsec].stackPrmryOff <= catalog[0].Nmeasure, "stackPrmryOff out of range"); - - meas = bestEntryMax; + off_t meas = bestEntryMax; myAssert (meas >= 0, "this should not happen"); @@ -687,44 +699,15 @@ // data for which the associated image has not been loaded (probably because of // overlaps). Msys + measure.Mcal is our best guess of the true magnitude - Mmos = Mgrid = 0; - McalPSF = measure[meas].McalPSF; // check that this is zero for loaded REF value + McalPSF = measure[meas].McalPSF; // check that this is zero for loaded REF value McalAPER = McalPSF; // check that this is zero for loaded REF value + Mmos = 0.0; + Mgrid = 0.0; } else { - McalPSF = getMcal (meas, cat, MAG_CLASS_PSF); - if (isnan(McalPSF)) SKIP_THIS_MEAS_STACK(Ncal); // XXX really skip? - + McalPSF = getMcal (meas, cat, MAG_CLASS_PSF); McalAPER = getMcal (meas, cat, MAG_CLASS_KRON); - if (isnan(McalAPER)) SKIP_THIS_MEAS_STACK(Ncal); // XXX really skip? - - Mmos = getMmos (meas, cat); - if (isnan(Mmos)) SKIP_THIS_MEAS_STACK(Nmos); - Mgrid = getMgrid (meas, cat); - if (isnan(Mgrid)) SKIP_THIS_MEAS_STACK(Ngrid); + Mmos = getMmos (meas, cat); + Mgrid = getMgrid (meas, cat); } - // NOTE: negative and insignificant flux values are allowed, but not NAN flux values - float Finst = PhotFluxInst (&measure[meas], MAG_CLASS_PSF); - if (isnan(Finst)) SKIP_THIS_MEAS_STACK(Ninst); - - // data quality assessment - isBad |= (measure[meas].photFlags & code->photomBadMask); - isBad |= (measure[meas].psfQF < 0.85); - isBad |= isnan(measure[meas].psfQF); - isBad |= measure[meas].dM > 0.2; // S/N < 5.0 - - isSuspect |= (measure[meas].photFlags & code->photomPoorMask); - isSuspect |= (measure[meas].psfQFperf < 0.85); - - if (!isSuspect && !isBad) { - NstackGood ++; - } - if (isSuspect && !isBad) { - NstackSuspect ++; - } - - // we are now populating stackDetectID not stack Image ID in secfilt - // ID = (stackPrimaryMeasureMin >= 0) ? stackPrimaryIDmin : stackCenterIDmin; - // ID = measure[meas].extID; // for the stack, this is the stackDetectID - // ** Here is the math to relate mag,zp to flux @@ -780,16 +763,17 @@ // Jy to AB mags - secfilt[Nsec].MpsfStk = (measure[meas].FluxPSF > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FpsfStk) : NAN; + secfilt[Nsec].MpsfStk = (measure[meas].FluxPSF > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FpsfStk) : NAN; secfilt[Nsec].MkronStk = (measure[meas].FluxKron > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FkronStk) : NAN; - secfilt[Nsec].MapStk = (measure[meas].FluxAp > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FapStk) : NAN; + secfilt[Nsec].MapStk = (measure[meas].FluxAp > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FapStk) : NAN; // record the measurement which gave the best value - secfilt[Nsec].stackBestOff = (bestEntryMax > -1) ? bestEntryMax + measureOffset; + secfilt[Nsec].stackBestOff = (bestEntryMax > -1) ? -1 : bestEntryMax + measureOffset; myAssert (secfilt[Nsec].stackBestOff <= catalog[0].Nmeasure, "stackBestOff out of range"); - secfilt[Nsec].Nstack = NstackMeas; - secfilt[Nsec].NstackDet = NstackDet; - - // this is the measurement used by secfilt[] + // record the selected primary measurement + secfilt[Nsec].stackPrmryOff = (primaryEntryMax == -1) ? -1 : primaryEntryMax + measureOffset; + myAssert (secfilt[Nsec].stackPrmryOff <= catalog[0].Nmeasure, "stackPrmryOff out of range"); + + // this is the selected measurement used by secfilt[] measure[meas].dbFlags |= ID_MEAS_STACK_PHOT_SRC; if (Nprimary) { @@ -801,7 +785,8 @@ // stack measurement used for 'best' was a detection (not forced from the other bands) - if ((measure[meas].photFlags2 & 0x00000004) == 0) { + if ((measure[bestEntryMax].photFlags2 & 0x00000004) == 0) { secfilt[Nsec].flags |= ID_SECF_STACK_BESTDET; } + // stack measurement used for 'primary' was a detection (not forced from the other bands) if ((primaryEntryMax >= 0) && ((measure[primaryEntryMax].photFlags2 & 0x00000004) == 0)) { @@ -810,4 +795,27 @@ secfilt[Nsec].flags |= ID_SECF_HAS_PS1_STACK; + + // NOTE: negative and insignificant flux values are allowed, but not NAN flux values + // float Finst = PhotFluxInst (&measure[meas], MAG_CLASS_PSF); + + // data quality assessment + isBad |= (measure[meas].photFlags & photcodes[Ns][0].photomBadMask); + isBad |= (measure[meas].psfQF < 0.85); + isBad |= isnan(measure[meas].psfQF); + isBad |= measure[meas].dM > 0.2; // S/N < 5.0 + + isSuspect |= (measure[meas].photFlags & photcodes[Ns][0].photomPoorMask); + isSuspect |= (measure[meas].psfQFperf < 0.85); + + if (!isSuspect && !isBad) { + NstackGood ++; + } + if (isSuspect && !isBad) { + NstackSuspect ++; + } + + secfilt[Nsec].Nstack = NstackMeas; + secfilt[Nsec].NstackDet = NstackDet; + } // Nsecfilt loop @@ -820,5 +828,5 @@ int PrimaryIsBest = TRUE; - for (Ns = 0; Ns < Nphotcodes; Ns++) { + for (int Ns = 0; Ns < Nphotcodes; Ns++) { if (!(secfilt[Ns].flags & ID_SECF_HAS_PS1_STACK)) continue; // no stack detection in PS1, nothing is best if (secfilt[Ns].flags & ID_SECF_STACK_PRIMARY) continue; // primary stack detection is best