Changeset 40248

Timestamp:

Dec 6, 2017, 8:25:46 AM (9 years ago)

Author:

eugene

Message:

trying to fix up logic for stack BEST & PRIMARY (review)

File:

• branches/eam_branches/ohana.20170822/src/relphot/src/setMrelCatalog.c (modified) (10 diffs)

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

@@ -540 +540 @@
   SecFilt *secfilt = &catalog[0].secfilt[ave*Nsecfilt];
 
-  off_t k;
-
   float McalPSF = 0, McalAPER = 0, Mmos = 0, Mgrid = 0;
 
@@ -558 +556 @@
   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 +565 @@
     int haveStack = FALSE;
 
-    float psfQFbest = 0.0;
-    float SNbest = 0.0;
-
-    off_t stackBestMeasure = -1;
     off_t stackPrimaryMeasure = -1;
 
@@ -577 +570 @@
     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 +583 @@
 
     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 +593 @@
       if (code->equiv != thisCode) { continue; }
 
-      Nstack ++;
+      NstackMeas ++;
       if ((measure[k].photFlags2 & 0x00000004) == 0) NstackDet ++;
 
@@ -603 +602 @@
       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 +722 @@
       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 +758 @@
     // 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 +768 @@
 
     // 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;
     }