Changeset 37113

Timestamp:

Jul 24, 2014, 11:14:51 AM (12 years ago)

Author:

eugene

Message:

stack primary bit does not depend on data quality; define stack best, stack primary bits correctly; clean_measures does not touch warp or stack; mean warp limited to same skycell as best stack

Location:

branches/eam_branches/ipp-20140717/Ohana/src/relphot

Files:

• include/relphot.h (modified) (1 diff)
• src/ImageOps.c (modified) (1 diff)
• src/StarOps.c (modified) (2 diffs)
• src/bcatalog.c (modified) (1 diff)
• src/setMrelCatalog.c (modified) (14 diffs)

branches/eam_branches/ipp-20140717/Ohana/src/relphot/include/relphot.h

@@ -440 +440 @@
 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);

branches/eam_branches/ipp-20140717/Ohana/src/relphot/src/ImageOps.c

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

branches/eam_branches/ipp-20140717/Ohana/src/relphot/src/StarOps.c

@@ -1086 +1086 @@
           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 +1137 @@
           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);

branches/eam_branches/ipp-20140717/Ohana/src/relphot/src/bcatalog.c

@@ -78 +78 @@
         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]);
 

branches/eam_branches/ipp-20140717/Ohana/src/relphot/src/setMrelCatalog.c

@@ -441 +441 @@
 
 // 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 +471 @@
     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 +481 @@
     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 +492 @@
       // 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 +503 @@
         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 +529 @@
       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 +558 @@
 
     // 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 +608 @@
     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 +646 @@
 
   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 +668 @@
 
   // 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 +684 @@
     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 +703 @@
       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 +733 @@
       // 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 +752 @@
       wpsflist[Npsf] = 1.0;
       Npsf ++;
+
+      measure[k].dbFlags |= ID_MEAS_WARP_USED;
 
       float Fap = PhotFluxCat (&measure[k], MAG_CLASS_APER);