Changeset 40249

Timestamp:

Dec 7, 2017, 5:54:38 AM (9 years ago)

Author:

eugene

Message:

rework of the assessment of stack detection quality

Location:

branches/eam_branches/ohana.20170822/src

Files:

• libdvo/include/dvo.h (modified) (1 diff)
• relphot/src/setMrelCatalog.c (modified) (15 diffs)

branches/eam_branches/ohana.20170822/src/libdvo/include/dvo.h

@@ -248 +248 @@
 
   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;
 

branches/eam_branches/ohana.20170822/src/relphot/src/setMrelCatalog.c

@@ -565 +565 @@
     int haveStack = FALSE;
 
-    off_t stackPrimaryMeasure = -1;
-
     int isBad = FALSE;
     int isSuspect = FALSE;
@@ -572 +570 @@
     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 +581 @@
     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 +608 @@
 
       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 +617 @@
       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 +624 @@
 
       // 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 +632 @@
       }
 
-      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 +644 @@
 
       // 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 +665 @@
 
       // 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 +687 @@
     // 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 +699 @@
       // 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 +763 @@
 
     // 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 +785 @@
 
     // 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 +795 @@
 
     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 +828 @@
 
   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