Index: /branches/eam_branches/ohana.20170822/src/relphot/src/setMrelCatalog.c =================================================================== --- /branches/eam_branches/ohana.20170822/src/relphot/src/setMrelCatalog.c (revision 40247) +++ /branches/eam_branches/ohana.20170822/src/relphot/src/setMrelCatalog.c (revision 40248) @@ -540,6 +540,4 @@ SecFilt *secfilt = &catalog[0].secfilt[ave*Nsecfilt]; - off_t k; - float McalPSF = 0, McalAPER = 0, Mmos = 0, Mgrid = 0; @@ -558,6 +556,5 @@ int haveStackObject = FALSE; - int Ns; - for (Ns = 0; Ns < Nphotcodes; Ns++) { + for (int Ns = 0; Ns < Nphotcodes; Ns++) { int thisCode = photcodes[Ns][0].code; @@ -568,8 +565,4 @@ int haveStack = FALSE; - float psfQFbest = 0.0; - float SNbest = 0.0; - - off_t stackBestMeasure = -1; off_t stackPrimaryMeasure = -1; @@ -577,6 +570,12 @@ int isSuspect = FALSE; - int Nstack = 0; // number for this photcode - int NstackDet = 0; // number for this photcode + int NstackMeas = 0; // number of stack measurements for this photcode + int NstackDet = 0; // number of stack detections for this photcode (not forced) + + int bestLevelMax = 0; + off_t bestEntryMax = -1; + float bestValueMax = 0.0; + + int Nprimary = 0; secfilt[Nsec].stackBestOff = -1; @@ -584,5 +583,5 @@ off_t meas = measureOffset; - for (k = 0; k < Nmeasure; k++, meas++) { + for (off_t k = 0; k < Nmeasure; k++, meas++) { // only examine gpc1 stack data @@ -594,5 +593,5 @@ if (code->equiv != thisCode) { continue; } - Nstack ++; + NstackMeas ++; if ((measure[k].photFlags2 & 0x00000004) == 0) NstackDet ++; @@ -603,74 +602,117 @@ haveStack = TRUE; haveStackObject = TRUE; + + XXX ??? if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad); + + int isPrimary = FALSE; - // ** 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 + // ** is this a PRIMARY measurement? (there may be more than one) + + // 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; - myAssert (secfilt[Nsec].stackPrmryOff <= catalog[0].Nmeasure, "stackPrmryOff out of range"); - } - - // ** now choose the BEST measurements (may also be PRIMARY) - - // ensure that we at least have a single best measure - if (stackBestMeasure == -1) stackBestMeasure = k; - - // choose the best psfQFperf value for the BEST measurement - if (isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > psfQFbest)) { - psfQFbest = measure[k].psfQFperf; - stackBestMeasure = k; - } - // ... UNLESS psfQFperf > 0.98 for the primary, in which case just use the primary. - if ((measure[k].dbFlags & ID_MEAS_STACK_PRIMARY) && isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > 0.98)) { - psfQFbest = 1000; // force this to be the best entry - stackBestMeasure = k; - } - - if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad); - + isPrimary = TRUE; + Nprimary ++; + } + + int psfQFperfAboveLimit = isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > 0.95); + + // ** determine the BEST level + int bestLevel = 0; + if (!isPrimary && !psfQFperfAboveLimit) bestLevel = 1; + if ( isPrimary && !psfQFperfAboveLimit) bestLevel = 2; + if (!isPrimary && psfQFperfAboveLimit) bestLevel = 3; + if ( isPrimary && psfQFperfAboveLimit) bestLevel = 4; + + // 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; + + // if we have not reached this level before, set the new level + if (bestLevelMax < bestLevel) { + bestValueMax = bestValue; + bestLevelMax = bestLevel; + bestEntryMax = k; + } + // if we have reached this level before, set the new value if we beat the old one + if ((bestLevelMax == bestLevel) && (bestValueMax < bestValue)) { + bestValueMax = bestValue; + bestEntryMax = k; + } + myAssert (bestEntryMax > -1, "this should not happen"); + + // ** determine the PRIMARY level + int primaryLevel = 0; + if ( isPrimary && !psfQFperfAboveLimit) primaryLevel = 1; + if ( isPrimary && psfQFperfAboveLimit) primaryLevel = 2; + + // 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; + if (primaryLevel == 1) { primaryValue = measure[k].psfQFperf; } + if (primaryLevel == 2) { primaryValue = 1.0 / measure[k].dM; } + + // if we have not reached this level before, set the new level + if (isPrimary && (primaryLevelMax < primaryLevel)) { + primaryValueMax = primaryValue; + primaryLevelMax = primaryLevel; + primaryEntryMax = k; + } + // if we have reached this level before, set the new value if we beat the old one + if (isPrimary && (primaryLevelMax == primaryLevel) && (primaryValueMax < primaryValue)) { + primaryValueMax = primaryValue; + primaryEntryMax = k; + } + } + + // if we do not have a stack measurement for this photcode, skip everything below + if (!haveStack) continue; + + // 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; + myAssert (meas >= 0, "this should not happen"); + + // measurements without an image are either external reference photometry or + // data for which the associated image has not been loaded (probably because of + // overlaps). we only want measurements associated with stack images in this loop + + // match measurement to its image (this is just a check, right?) + if (getImageEntry (meas, cat) < 0) { // measurements without an image are either external reference photometry or // data for which the associated image has not been loaded (probably because of - // overlaps). we only want measurements associated with stack images in this loop - - // match measurement to its image (this is just a check, right?) - if (getImageEntry (meas, cat) < 0) { - // measurements without an image are either external reference photometry or - // 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[k].McalPSF; // check that this is zero for loaded REF value - McalAPER = McalPSF; // check that this is zero for loaded REF value - } else { - McalPSF = getMcal (meas, cat, MAG_CLASS_PSF); - if (isnan(McalPSF)) SKIP_THIS_MEAS_STACK(Ncal); // XXX really skip? - - 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); - } - - // NOTE: negative and insignificant flux values are allowed, but not NAN flux values - float Finst = PhotFluxInst (&measure[k], MAG_CLASS_PSF); - if (isnan(Finst)) SKIP_THIS_MEAS_STACK(Ninst); - - // data quality assessment - isBad |= (measure[k].photFlags & code->photomBadMask); - isBad |= (measure[k].psfQF < 0.85); - isBad |= isnan(measure[k].psfQF); - isBad |= measure[k].dM > 0.2; // S/N < 5.0 - - isSuspect |= (measure[k].photFlags & code->photomPoorMask); - isSuspect |= (measure[k].psfQFperf < 0.85); - } - - if (!haveStack) continue; + // 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 + McalAPER = McalPSF; // check that this is zero for loaded REF value + } else { + McalPSF = getMcal (meas, cat, MAG_CLASS_PSF); + if (isnan(McalPSF)) SKIP_THIS_MEAS_STACK(Ncal); // XXX really skip? + + 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); + } + + // 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) { @@ -680,17 +722,10 @@ NstackSuspect ++; } - - // measurements which are bad will not have a valid stack entry and are skipped - k = (stackBestMeasure >= 0) ? stackBestMeasure : stackPrimaryMeasure; - if (k < 0) continue; - + // we are now populating stackDetectID not stack Image ID in secfilt // ID = (stackPrimaryMeasureMin >= 0) ? stackPrimaryIDmin : stackCenterIDmin; - // ID = measure[k].extID; // for the stack, this is the stackDetectID - - // get the zero point for the selected image - // Use a different zero point for the PSF-like and APERTURE-like magnitudes - float zpPSF = PhotZeroPoint (&measure[k], &average[0], &secfilt[0]) - (McalPSF + Mmos + Mgrid); - float zpAPER = PhotZeroPoint (&measure[k], &average[0], &secfilt[0]) - (McalAPER + Mmos + Mgrid); + // ID = measure[meas].extID; // for the stack, this is the stackDetectID + + // ** Here is the math to relate mag,zp to flux // flux_cgs : erg sec^1 cm^-2 Hz^-1 @@ -723,4 +758,9 @@ // flux_Jy = flux_inst * ten(-0.4*ZP + 3.56) + // get the zero point for the selected image + // Use a different zero point for the PSF-like and APERTURE-like magnitudes + float zpPSF = PhotZeroPoint (&measure[meas], &average[0], &secfilt[0]) - (McalPSF + Mmos + Mgrid); + float zpAPER = PhotZeroPoint (&measure[meas], &average[0], &secfilt[0]) - (McalAPER + Mmos + Mgrid); + // zpFactor to go from instrumental flux to Janskies float zpFactorPSF = pow(10.0, -0.4*zpPSF + 3.56); @@ -728,38 +768,42 @@ // need to put in AB mag factor to get to Janskies (or uJy?) - secfilt[Nsec].FpsfStk = zpFactorPSF * measure[k].FluxPSF; - secfilt[Nsec].dFpsfStk = zpFactorPSF * measure[k].dFluxPSF; - secfilt[Nsec].FkronStk = zpFactorAPER * measure[k].FluxKron; - secfilt[Nsec].dFkronStk = zpFactorAPER * measure[k].dFluxKron; - secfilt[Nsec].FapStk = zpFactorAPER * measure[k].FluxAp; + secfilt[Nsec].FpsfStk = zpFactorPSF * measure[meas].FluxPSF; + secfilt[Nsec].dFpsfStk = zpFactorPSF * measure[meas].dFluxPSF; + secfilt[Nsec].FkronStk = zpFactorAPER * measure[meas].FluxKron; + secfilt[Nsec].dFkronStk = zpFactorAPER * measure[meas].dFluxKron; + secfilt[Nsec].FapStk = zpFactorAPER * measure[meas].FluxAp; // NOTE: for PV3, apFluxErr is broken (see pmSourcePhotometry.c:245). we are going to // use PSF flux error instead here: - // secfilt[Nsec].dFapStk = zpFactorAPER * measure[k].dFluxAp; - secfilt[Nsec].dFapStk = zpFactorAPER * measure[k].dFluxPSF; + // secfilt[Nsec].dFapStk = zpFactorAPER * measure[meas].dFluxAp; + secfilt[Nsec].dFapStk = zpFactorAPER * measure[meas].dFluxPSF; // 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].stackBestOff = k + measureOffset; + 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; + + // record the measurement which gave the best value + secfilt[Nsec].stackBestOff = (bestEntryMax > -1) ? bestEntryMax + measureOffset; myAssert (secfilt[Nsec].stackBestOff <= catalog[0].Nmeasure, "stackBestOff out of range"); - secfilt[Nsec].Nstack = Nstack; + secfilt[Nsec].Nstack = NstackMeas; secfilt[Nsec].NstackDet = NstackDet; // this is the measurement used by secfilt[] - measure[k].dbFlags |= ID_MEAS_STACK_PHOT_SRC; - if (k == stackPrimaryMeasure) { + measure[meas].dbFlags |= ID_MEAS_STACK_PHOT_SRC; + if (Nprimary) { secfilt[Nsec].flags |= ID_SECF_STACK_PRIMARY; + if (Nprimary > 1) { + secfilt[Nsec].flags |= ID_SECF_STACK_PRIMARY_MULTIPLE; + } } // stack measurement used for 'best' was a detection (not forced from the other bands) - if ((measure[k].photFlags2 & 0x00000004) == 0) { + if ((measure[meas].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 ((stackPrimaryMeasure >= 0) && ((measure[stackPrimaryMeasure].photFlags2 & 0x00000004) == 0)) { + if ((primaryEntryMax >= 0) && ((measure[primaryEntryMax].photFlags2 & 0x00000004) == 0)) { secfilt[Nsec].flags |= ID_SECF_STACK_PRIMDET; }