Changeset 34772 for branches/czw_branch/20120906/Ohana/src/relphot/src/StarOps.c

Timestamp:

Dec 6, 2012, 6:32:19 PM (14 years ago)

Author:

watersc1

Message:

merge from trunk. Preliminary versions of stacktool/warptool updates and regenerate_background.pl script. The warp code is finished and tested, and I still need to get the stack version resolved.

Location:

branches/czw_branch/20120906/Ohana

Files:

• . (modified) (1 prop)
• src/relphot/src/StarOps.c (modified) (26 diffs)

branches/czw_branch/20120906/Ohana/src/relphot/src/StarOps.c

@@ -11 +11 @@
   int Nmos;
   int Ngrid;
-  double *list;
-  double *dlist;
-  double *wlist;
-  double *aplist;
-  double *kronlist;
-  double *dkronlist;
+
+  // NOTE: the following arrays are (possibly) pre-allocated and carried down to each
+  // thread.  The first 3 (list, dlist, wlist) are used in all relphot analyses; the
+  // others are only used on the final output steps.
+
+  double *list;               // list of measure.mag values for a given star
+  double *dlist;              // mag errors for a star
+  double *wlist;              // weights to use for mean mags
+  double *aplist;             // ap mags for a star
+  double *kronlist;           // kron mags for a star
+  double *dkronlist;          // kron mag errors 
+
+  double *psfqf_list;         // psfqf for all filters
+  double *psfqfperf_list;     // psfqfperf for all filters
+  double *stargal_list;       // stargal for all filters
 } SetMrelInfo;
 
@@ -168 +177 @@
   ALLOCATE (results.dkronlist, double, Nmax);
 
+  ALLOCATE (results.psfqf_list, double, Nmax);
+  ALLOCATE (results.psfqfperf_list, double, Nmax);
+  ALLOCATE (results.stargal_list, double, Nmax);
+
   for (i = 0; i < Ncatalog; i++) {
     setMrel_catalog  (catalog, i, pass, flatcorr, &results, Nsecfilt); // XXX add arguments as needed for options
@@ -178 +191 @@
   free (results.kronlist);
   free (results.dkronlist);
+
+  free (results.psfqf_list);
+  free (results.psfqfperf_list);
+  free (results.stargal_list);
   return (TRUE);
 }
@@ -315 +332 @@
   liststats_setmode (&kronstats, STATMODE);
 
+  // StatType psfstats;
+  // liststats_setmode (&psfstats, "MEAN");
+
   double *list      = results->list;
   double *dlist     = results->dlist;
@@ -322 +342 @@
   double *dkronlist = results->dkronlist;
 
+  double *psfqf_list     = results->psfqf_list;
+  double *psfqfperf_list = results->psfqfperf_list;
+  double *stargal_list   = results->stargal_list;
+
   SetMrelInfoInit (results, FALSE); // do not allocate list,dlist,wlist arrays
 
@@ -327 +351 @@
 
   char *primaryCell = NULL;
-  ALLOCATE (primaryCell, char, DVO_MAX_PATH);
+  if (isSetMrelFinal) {
+    ALLOCATE (primaryCell, char, DVO_MAX_PATH);
+  }
 
   for (j = 0; j < catalog[Nc].Naverage; j++) {
-    // XXX accumulate all secfilt values in a single pass?
-
     // option for a test print
     if (FALSE && (catalog[Nc].average[j].objID == 0x7146) && (catalog[Nc].average[j].catID == 0x49d8)) {
@@ -338 +362 @@
     }
 
-    BoundaryTreePrimaryCell(primaryCell, catalog[Nc].average[j].R, catalog[Nc].average[j].D);
+    if (isSetMrelFinal) {
+      BoundaryTreePrimaryCell(primaryCell, catalog[Nc].average[j].R, catalog[Nc].average[j].D);
+    }
 
     int GoodPS1 = FALSE;
@@ -346 +372 @@
     int NextPS1 = 0;
     int NpsfPS1 = 0;
+
+    int Nmeas = 0; // number of meas used for psfqf, psfqfperf, and stargal
+
+    int NstackGood = 0;
+    int NstackSuspect = 0;
+
+    int haveStackObject = FALSE;
 
     int Ns;
@@ -365 +398 @@
       int haveSynth = FALSE;
       int haveStack = FALSE;
+      int havePS1   = FALSE;
 
       // need to find the measurement closest to the center of its skycell, as well as the
@@ -370 +404 @@
 
       float stackCenterOffsetMin = 1e9;
-      int stackCenterIDmin = -1;
       off_t stackCenterMeasureMin = -1;
 
       float stackPrimaryOffsetMin = 1e9;
-      int stackPrimaryIDmin = -1;
       off_t stackPrimaryMeasureMin = -1;
 
@@ -384 +416 @@
       int minUbercalDist = 1000;
     
+      int isBad = FALSE;
+      int isSuspect = FALSE;
+
       N = 0;
       m = catalog[Nc].averageT[j].measureOffset;
@@ -393 +428 @@
         if (code->equiv != thisCode) { continue; }
         Ncode ++;
+
+        catalog[Nc].measure[m].dbFlags &= ~ID_MEAS_STACK_PRIMARY;
 
         if (catalog[Nc].measureT[m].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); 
@@ -414 +451 @@
         if (isnan(Msys)) SKIP_THIS_MEAS(Nsys);
 
-        list[N] = Msys - Mcal - Mmos - Mgrid;
-
         int myUbercalDist = getUbercalDist(m,Nc);
         minUbercalDist = MIN(minUbercalDist, myUbercalDist);
 
         if (isSetMrelFinal) {
-          float Map = PhotAper (&catalog[Nc].measure[m]);
-          aplist[N] = Map - Mcal - Mmos - Mgrid;
-
-          float Mkron = PhotKron (&catalog[Nc].measure[m]);
-          kronlist[N] = Mkron - Mcal - Mmos - Mgrid;
-          dkronlist[N] = catalog[Nc].measure[m].dMkron;
-
           // special options for PS1 data
           if ((catalog[Nc].measure[m].photcode >= 10000) && (catalog[Nc].measure[m].photcode <= 10500)) {
             // count the extended detections
             if (!isnan(catalog[Nc].measure[m].Map)) {
-              if (catalog[Nc].measure[m].M - catalog[Nc].measure[m].Map > 0.5) {
+              float dMagAp = catalog[Nc].measure[m].M - catalog[Nc].measure[m].Map;
+              float SigmaAp = hypot(0.1, 2.5*catalog[Nc].measure[m].dM);
+              // XXX this is still quite ad hoc, but at least it:
+              // (a) converges to 0.1 mag offset at the bright end
+              // (b) converges to 0.5 mag offset at the faint end (dM = 0.2)
+              if (dMagAp > SigmaAp) {
                 Next ++;
                 NextPS1 ++;
@@ -438 +471 @@
               }
             }
+            havePS1 = TRUE;
           }
+
           // gpc1 stack data
           if ((catalog[Nc].measure[m].photcode >= 11000) && (catalog[Nc].measure[m].photcode <= 11400)) {
-            // if (pass < 2) continue;
             haveStack = TRUE;
+            haveStackObject = TRUE;
 
             unsigned int stackImageID;
@@ -452 +487 @@
               if (stackPrimaryOffset < stackPrimaryOffsetMin) {
                 stackPrimaryOffsetMin = stackPrimaryOffset;
-                stackPrimaryIDmin = stackImageID;
+                // stackPrimaryIDmin = stackImageID;
                 stackPrimaryMeasureMin = m;
               }
@@ -461 +496 @@
             if (stackCenterOffset < stackCenterOffsetMin) {
               stackCenterOffsetMin = stackCenterOffset;
-              stackCenterIDmin = stackImageID;
+              // stackCenterIDmin = stackImageID;
               stackCenterMeasureMin = m;
             }
+
+            // data quality assessment
+            isBad |= (catalog[Nc].measure[m].photFlags & code->photomBadMask);
+            isBad |= (catalog[Nc].measure[m].psfQF < 0.85);
+            isBad |= isnan(catalog[Nc].measure[m].psfQF);
+            isBad |= catalog[Nc].measure[m].dM > 0.2; // S/N < 5.0
+
+            isSuspect |= (catalog[Nc].measure[m].photFlags & code->photomPoorMask);
+            isSuspect |= (catalog[Nc].measure[m].psfQFperf < 0.85);
+
+            // mean psfQF, psfQFperf, and stargal values are calculated per object (all filters, stacks and exposures)
+            psfqf_list[Nmeas]     = catalog[Nc].measure[m].psfQF;
+            psfqfperf_list[Nmeas] = catalog[Nc].measure[m].psfQFperf;
+            stargal_list[Nmeas]   = catalog[Nc].measure[m].extNsigma;
+            Nmeas++;
+
+            continue; // do not include stack measurements in the per-exposure statistics
           }
 
@@ -478 +530 @@
                 // detections without one of these bits should only be used in PASS_1
                 SKIP_THIS_MEAS(Nbad);
-                continue;
               }
             }
@@ -512 +563 @@
               if (pass < 4) {
                 SKIP_THIS_MEAS(Nbad);
-                continue;
               }
               haveSynth = TRUE;
             }
           }
+
+          // Map (aplist) and Mkron (kronlist,dkronlist) are used to calculate mean mags per filter
+          float Map = PhotAper (&catalog[Nc].measure[m]);
+          aplist[N] = Map - Mcal - Mmos - Mgrid;
+
+          float Mkron = PhotKron (&catalog[Nc].measure[m]);
+          kronlist[N] = Mkron - Mcal - Mmos - Mgrid;
+          dkronlist[N] = catalog[Nc].measure[m].dMkron;
+
+          // mean psfQF, psfQFperf, and stargal values are calculated per object (all filters, stacks and exposures)
+          psfqf_list[Nmeas]     = catalog[Nc].measure[m].psfQF;
+          psfqfperf_list[Nmeas] = catalog[Nc].measure[m].psfQFperf;
+          stargal_list[Nmeas]   = catalog[Nc].measure[m].extNsigma;
+          Nmeas++;
         }
 
@@ -527 +591 @@
         dlist[N] = MAX (hypot(catalog[Nc].measureT[m].dM, code->photomErrSys), MIN_ERROR);
         wlist[N] = 1.0;
+        list[N]  = Msys - Mcal - Mmos - Mgrid;
 
         // up-weight the ubercal values (or convergence can take a long time...)
@@ -549 +614 @@
       if (haveStack && (N > 1)) {
         // fprintf (stderr, "multiple stack values for %10.6f %10.6f\n", catalog[Nc].averageT[j].R, catalog[Nc].averageT[j].D);
+      }
+
+      if (haveStack) {
+        if (!isSuspect && !isBad) {
+          NstackGood ++;
+        }
+        if (isSuspect && !isBad) {
+          NstackSuspect ++;
+        }
       }
 
@@ -604 +678 @@
         if (haveStack) {
           m  = (stackPrimaryMeasureMin >= 0) ? stackPrimaryMeasureMin : stackCenterMeasureMin;
-          ID = (stackPrimaryMeasureMin >= 0) ? stackPrimaryIDmin      : stackCenterIDmin;
+
+          // we are now populating stackDetectID not stack Image ID in secfilt
+          // ID = (stackPrimaryMeasureMin >= 0) ? stackPrimaryIDmin      : stackCenterIDmin;
+          ID = catalog[Nc].measure[m].extID; // for the stack, this is the stackDetectID
 
           // get the zero point for the selected image
-          float zp = Mcal + Mmos + Mgrid + PhotZeroPoint (&catalog[Nc].measure[m], &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt]);
+          float zp = PhotZeroPoint (&catalog[Nc].measure[m], &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt]) - (Mcal + Mmos + Mgrid);
 
           // flux_cgs : erg sec^1 cm^-2 Hz^-1
@@ -647 +724 @@
           catalog[Nc].secfilt[Nsecfilt*j+Nsec].dFluxKron = zpFactor * catalog[Nc].measure[m].dFluxKron;
 
-          catalog[Nc].secfilt[Nsecfilt*j+Nsec].stackID   = ID;
+          catalog[Nc].secfilt[Nsecfilt*j+Nsec].stackDetectID = ID;
+
+          catalog[Nc].measure[m].dbFlags |= ID_MEAS_STACK_PHOT_SRC;
+          if (stackPrimaryMeasureMin >= 0) {
+            catalog[Nc].measure[m].dbFlags |= ID_MEAS_STACK_PRIMARY;
+          }
         }
 
@@ -658 +740 @@
           case 0:
             catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_0;
-            GoodPS1 = TRUE;
+            if (havePS1) GoodPS1 = TRUE;
             break;
           case 1:
             catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_1;
-            GoodPS1 = TRUE;
+            if (havePS1) GoodPS1 = TRUE;
             break;
           case 2:
             catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_2;
-            GoodPS1 = TRUE;
+            if (havePS1) GoodPS1 = TRUE;
             break;
           case 3:
@@ -678 +760 @@
           catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_SYNTH;
         }       
+        if (haveStack) {
+          catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_HAS_STACK;
+        }       
+        if (havePS1) {
+          catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_HAS_PS1;
+        }       
         if (haveUbercal) {
           catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_UBERCAL;
@@ -685 +773 @@
 
     if (isSetMrelFinal) {
-      DVOAverageFlags flagBits = ID_OBJ_EXT | ID_OBJ_EXT_ALT | ID_OBJ_GOOD | ID_OBJ_GOOD_ALT;
-
-      // we attempt to set a few flags here; reset those bits before trying:
-      catalog[Nc].average[j].flags &= ~flagBits;
+      if (pass == 0) {
+        DVOAverageFlags flagBits = ID_OBJ_EXT | ID_OBJ_EXT_ALT | ID_OBJ_GOOD | ID_OBJ_GOOD_ALT;
+        // we attempt to set a few flags here; reset those bits before trying:
+        catalog[Nc].average[j].flags &= ~flagBits;
+      }
 
       if (NextPS1 && (NextPS1 > NpsfPS1)) {
@@ -701 +790 @@
       if (Good2MASS) {
         catalog[Nc].average[j].flags |= ID_OBJ_GOOD_ALT;
+      }
+
+      if (NstackGood >= 2) {
+        catalog[Nc].average[j].flags |= ID_OBJ_GOOD_STACK;
+      } else if (NstackGood + NstackSuspect >= 2) {
+        catalog[Nc].average[j].flags |= ID_OBJ_SUSPECT_STACK;
+      } else if (haveStackObject) {
+        catalog[Nc].average[j].flags |= ID_OBJ_BAD_STACK;
+      }
+
+      // use NAN values for catalog[Nc].average[j].psfQF as signal
+      if (!isfinite(catalog[Nc].average[j].psfQF)) {
+        // do something here
+        float psfQFmax     = 0.0;
+        float psfQFperfmax = 0.0;
+        float stargalmax   = 0.0;
+        int seq;
+        for (seq = 0; seq < Nmeas; seq++) {
+          if (isfinite(psfqf_list[seq])) {
+            psfQFmax     = MAX (psfqf_list[seq], psfQFmax);
+          }
+          if (isfinite(psfqfperf_list[seq])) {
+            psfQFperfmax = MAX (psfqfperf_list[seq], psfQFperfmax);
+          }
+          if (isfinite(stargal_list[seq])) {
+            stargalmax   = MAX (stargal_list[seq], stargalmax);
+          }
+        }
+
+        catalog[Nc].average[j].psfQF     = psfQFmax;
+        catalog[Nc].average[j].psfQFperf = psfQFperfmax;
+        catalog[Nc].average[j].stargal   = stargalmax;
       }
     }