Changeset 34260 for trunk/Ohana/src/relphot

Timestamp:

Jul 31, 2012, 4:02:00 PM (14 years ago)

Author:

eugene

Message:

add PS1_V4 schema; add set mean fluxes based on stacks; define the primary image ID (based on distance from center or boundary tree) and set mean flux from stacks based on that info; fix (uniquify) psps ID for stack detections; mean fluxes are in Jy; detection fluxes are instrumental counts per second

Location:

trunk/Ohana/src/relphot

Files:

• Makefile (modified) (2 diffs)
• include/relphot.h (modified) (3 diffs)
• src/BoundaryTreeOps.c (copied) (copied from branches/eam_branches/ipp-20120627/Ohana/src/relphot/src/BoundaryTreeOps.c )
• src/ImageOps.c (modified) (1 diff)
• src/StarOps.c (modified) (10 diffs)
• src/args.c (modified) (1 diff)

trunk/Ohana/src/relphot/Makefile

@@ -49 +49 @@
 $(SRC)/setExclusions.$(ARCH).o   \
 $(SRC)/setMrelFinal.$(ARCH).o    \
+$(SRC)/BoundaryTreeOps.$(ARCH).o         \
 $(SRC)/write_coords.$(ARCH).o
 
@@ -77 +78 @@
 $(SRC)/setExclusions.$(ARCH).o   \
 $(SRC)/setMrelFinal.$(ARCH).o    \
+$(SRC)/BoundaryTreeOps.$(ARCH).o         \
 $(SRC)/write_coords.$(ARCH).o
 

trunk/Ohana/src/relphot/include/relphot.h

@@ -111 +111 @@
 char        *BCATALOG;
 ModeType     MODE;
+
+char        *BOUNDARY_TREE;
 
 double MAG_LIM;
@@ -242 +244 @@
 float         getMrel             PROTO((Catalog *catalog, off_t meas, int cat));
 short         getUbercalDist      PROTO((off_t meas, int cat));
+float         getCenterOffset     PROTO((off_t meas, int cat, Measure *measure, unsigned int *myID));
 Image        *getimage            PROTO((off_t N));
 Image        *getimages           PROTO((off_t *N, off_t **LineNumber));
@@ -340 +343 @@
 int client_logger_message (char *format,...);
 
+int MatchImageName (off_t meas, int cat, char *name);
+
+int load_tree (char *treefile);
+int BoundaryTreePrimaryCell (char *primaryCellName, double ra, double dec);

trunk/Ohana/src/relphot/src/ImageOps.c

@@ -365 +365 @@
   distance = image[i].ubercalDist; // was dummy3 in structure
   return (distance);
+}
+
+// returns image.Mcal - ff(x,y)
+float getCenterOffset (off_t meas, int cat, Measure *measure, unsigned int *myID) {
+
+  off_t i;
+  float distance;
+
+  if (!MeasureToImage) return -1;
+
+  i = MeasureToImage[cat][meas];
+  if (i == -1) return (1000);
+
+  float Xcenter = 0.5*image[i].NX;
+  float Ycenter = 0.5*image[i].NY;
+
+  *myID = image[i].imageID;
+
+  distance = hypot (measure[0].Xccd - Xcenter, measure[0].Yccd - Ycenter);
+  return (distance);
+}
+
+// returns image.Mcal - ff(x,y)
+int MatchImageName (off_t meas, int cat, char *name) {
+
+  off_t i;
+
+  if (!name) return FALSE;
+  if (!name[0]) return FALSE;
+
+  if (!MeasureToImage) return FALSE;
+
+  i = MeasureToImage[cat][meas];
+  if (i == -1) return FALSE;
+
+  // this is a bit crude: stack image names are of the form 
+  // RINGS.V3.skycell.1495.027.sky.191211.stk.988232.cmf 
+  // the primaryCell has a name of the form RINGS.V3.skycell.1495
+
+  if (!strncmp(image[i].name, name, strlen(name))) return TRUE;
+  return FALSE;
 }
 

trunk/Ohana/src/relphot/src/StarOps.c

@@ -15 +15 @@
   double *wlist;
   double *aplist;
-  double *daplist;
+  double *kronlist;
+  double *dkronlist;
 } SetMrelInfo;
 
@@ -164 +165 @@
   SetMrelInfoInit (&results, TRUE); // allocates results->list,dlist,wlist
   ALLOCATE (results.aplist, double, Nmax);
-  ALLOCATE (results.daplist, double, Nmax);
+  ALLOCATE (results.kronlist, double, Nmax);
+  ALLOCATE (results.dkronlist, double, Nmax);
 
   for (i = 0; i < Ncatalog; i++) {
@@ -174 +176 @@
   SetMrelInfoFree (&results);
   free (results.aplist);
-  free (results.daplist);
+  free (results.kronlist);
+  free (results.dkronlist);
   return (TRUE);
 }
@@ -303 +306 @@
 int setMrel_catalog (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
 
-  off_t j, k, m;
+  off_t j, k, m, ID;
   int N;
   float Msys, Mcal, Mmos, Mgrid;
 
-  StatType stats, apstats;
+  StatType stats, apstats, kronstats;
   liststats_setmode (&stats, STATMODE);
   liststats_setmode (&apstats, STATMODE);
-
-  double *list    = results->list;
-  double *dlist   = results->dlist;
-  double *wlist   = results->wlist;
-  double *aplist  = results->aplist;
-  double *daplist = results->daplist;
+  liststats_setmode (&kronstats, STATMODE);
+
+  double *list      = results->list;
+  double *dlist     = results->dlist;
+  double *wlist     = results->wlist;
+  double *aplist    = results->aplist;
+  double *kronlist  = results->kronlist;
+  double *dkronlist = results->dkronlist;
 
   SetMrelInfoInit (results, FALSE); // do not allocate list,dlist,wlist arrays
 
   int isSetMrelFinal = (pass >= 0);
+
+  char *primaryCell = NULL;
+  ALLOCATE (primaryCell, char, DVO_MAX_PATH);
 
   for (j = 0; j < catalog[Nc].Naverage; j++) {
@@ -325 +333 @@
 
     // option for a test print
-    if (FALSE && (catalog[Nc].average[j].objID == 0x46a4) && (catalog[Nc].average[j].catID == 0xf40e)) {
+    if (FALSE && (catalog[Nc].average[j].objID == 0x7146) && (catalog[Nc].average[j].catID == 0x49d8)) {
       fprintf (stderr, "test obj\n");
       print_measure_set (&catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], catalog[Nc].measure);
     }
+
+    BoundaryTreePrimaryCell(primaryCell, catalog[Nc].average[j].R, catalog[Nc].average[j].D);
 
     int GoodPS1 = FALSE;
@@ -355 +365 @@
       int haveSynth = FALSE;
       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;
+      int stackCenterIDmin = -1;
+      off_t stackCenterMeasureMin = -1;
+
+      float stackPrimaryOffsetMin = 1e9;
+      int stackPrimaryIDmin = -1;
+      off_t stackPrimaryMeasureMin = -1;
 
       int forceSynth = FALSE;
@@ -401 +422 @@
           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
@@ -416 +441 @@
           // gpc1 stack data
           if ((catalog[Nc].measure[m].photcode >= 11000) && (catalog[Nc].measure[m].photcode <= 11400)) {
-            if (pass < 2) continue;
+            // if (pass < 2) continue;
             haveStack = TRUE;
+
+            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 (MatchImageName (m, Nc, primaryCell)) {
+              float stackPrimaryOffset = getCenterOffset (m, Nc, &catalog[Nc].measure[m], &stackImageID);
+              if (stackPrimaryOffset < stackPrimaryOffsetMin) {
+                stackPrimaryOffsetMin = stackPrimaryOffset;
+                stackPrimaryIDmin = stackImageID;
+                stackPrimaryMeasureMin = m;
+              }
+            }
+
+            // get the center distance for the generic case:
+            float stackCenterOffset = getCenterOffset (m, Nc, &catalog[Nc].measure[m], &stackImageID);
+            if (stackCenterOffset < stackCenterOffsetMin) {
+              stackCenterOffsetMin = stackCenterOffset;
+              stackCenterIDmin = stackImageID;
+              stackCenterMeasureMin = m;
+            }
           }
 
@@ -549 +595 @@
 
         // NOTE : use the modified weight for apmags as well as psf mags
-        liststats (aplist, daplist, wlist, N, &apstats);
-
+        liststats (aplist, dlist, wlist, N, &apstats);
         catalog[Nc].secfilt[Nsecfilt*j+Nsec].Map  = apstats.mean; 
+
+        liststats (kronlist, dkronlist, wlist, N, &kronstats);
+        catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mkron  = kronstats.mean; 
+        catalog[Nc].secfilt[Nsecfilt*j+Nsec].dMkron = kronstats.error; 
+
+        if (haveStack) {
+          m  = (stackPrimaryMeasureMin >= 0) ? stackPrimaryMeasureMin : stackCenterMeasureMin;
+          ID = (stackPrimaryMeasureMin >= 0) ? stackPrimaryIDmin      : stackCenterIDmin;
+
+          // 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]);
+
+          // flux_cgs : erg sec^1 cm^-2 Hz^-1
+          // mag_inst : -2.5 log (cts/sec)
+          // mag_inst : -2.5 log (flux_inst)
+          // flux_inst = ten(-0.4*mag_inst)
+
+          // mag_AB = -2.5 log (flux_cgs) - 48.6 (~by definition) [~Vega flux in V-band]
+          // flux_cgs = ten(-0.4*(mag_AB + 48.6))
+
+          // flux_AB : ten(-0.4*mag_AB)
+
+          // flux_cgs = ten(-0.4*48.6) * flux_AB
+          // flux_AB  = ten(+0.4*48.6) * flux_cgs
+            
+          // flux_Jy : flux_cgs * 10^23
+
+          // flux_AB = ten(+0.4*48.6) * ten(-23) * flux_Jy
+
+          // mag_AB = mag_inst + ZP
+
+          // flux_inst = ten(-0.4*(mag_AB - ZP)) = ten(0.4*ZP) * flux_AB
+
+          // flux_AB = flux_inst * ten(-0.4*ZP)
+
+          // flux_inst * ten(-0.4*ZP) = ten(+0.4*48.6 - 23) * flux_Jy
+
+          // flux_inst = flux_Jy * ten(0.4*ZP + 0.4*48.6 - 23)
+          // flux_inst = flux_Jy * ten(0.4*ZP - 3.56)
+          // flux_Jy = flux_inst * ten(-0.4*ZP + 3.56)
+
+          // zpFactor to go from instrumental flux to Janskies
+          float zpFactor = pow(10.0, -0.4*zp + 3.56);
+
+          // need to put in AB mag factor to get to Janskies (or uJy?)
+          catalog[Nc].secfilt[Nsecfilt*j+Nsec].FluxPSF   = zpFactor * catalog[Nc].measure[m].FluxPSF;  
+          catalog[Nc].secfilt[Nsecfilt*j+Nsec].dFluxPSF  = zpFactor * catalog[Nc].measure[m].dFluxPSF; 
+          catalog[Nc].secfilt[Nsecfilt*j+Nsec].FluxKron  = zpFactor * catalog[Nc].measure[m].FluxKron; 
+          catalog[Nc].secfilt[Nsecfilt*j+Nsec].dFluxKron = zpFactor * catalog[Nc].measure[m].dFluxKron;
+
+          catalog[Nc].secfilt[Nsecfilt*j+Nsec].stackID   = ID;
+        }
 
         // NOTE: for 2MASS measurements, Next should be 1, as should N
@@ -607 +704 @@
     }
   }
+  if (primaryCell) free (primaryCell);
   return (TRUE);
 }

trunk/Ohana/src/relphot/src/args.c

@@ -197 +197 @@
   }
 
+  if ((N = get_argument (argc, argv, "-boundary-tree"))) {
+    remove_argument (N, &argc, argv);
+    load_tree (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
   SHOW_PARAMS = FALSE;
   if ((N = get_argument (argc, argv, "-params"))) {