Changeset 42821 for trunk/Ohana/src/relphot

Timestamp:

May 8, 2025, 4:29:52 PM (14 months ago)

Author:

eugene

Message:

merge changes from eam_branches/ipp-20230313

Location:

trunk/Ohana

Files:

• . (modified) (1 prop)
• src/relphot (modified) (1 prop)
• src/relphot/include/relphot.h (modified) (3 diffs)
• src/relphot/src/ConfigInit.c (modified) (1 diff)
• src/relphot/src/ImageOps.c (modified) (3 diffs)
• src/relphot/src/StarOps.c (modified) (14 diffs)
• src/relphot/src/args.c (modified) (4 diffs)
• src/relphot/src/bcatalog.c (modified) (12 diffs)
• src/relphot/src/extra.c (modified) (3 diffs)
• src/relphot/src/fit1d_irls.c (modified) (2 diffs)
• src/relphot/src/initialize.c (modified) (1 diff)
• src/relphot/src/load_catalogs.c (modified) (2 diffs)
• src/relphot/src/reload_catalogs.c (modified) (1 diff)
• src/relphot/src/relphot.c (modified) (1 diff)

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

@@ -417 +417 @@
 int IS_DIFF_DB;
 
+double MAG_MIN;
 double MAG_LIM;
 double SIGMA_LIM;
+
+double STAR_MAX_PSF_KRON;
+double PSF_QF_MIN;
+
 double IMAGE_SCATTER;
 double IMAGE_OFFSET;
@@ -471 +476 @@
 int    RESET_FLATCORR;
 int    REPAIR_WARPS;
+int    USE_REF_EQUIV;
 int    PRESERVE_PS1;
 int    REQUIRE_PSFFIT;
@@ -486 +492 @@
 double IMAGE_GOOD_FRACTION;
 int    CALIBRATE_STACKS_AND_WARPS;
+int    CALIBRATE_WARPS_LIKE_CHIPS;
+int    CALIBRATE_STACKS_LIKE_CHIPS;
 
 int    KEEP_UBERCAL;

trunk/Ohana/src/relphot/src/ConfigInit.c

@@ -23 +23 @@
   if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file);
 
-  GetConfig (config, "MAG_LIM",                "%lf", 0, &MAG_LIM);
-  GetConfig (config, "SIGMA_LIM",              "%lf", 0, &SIGMA_LIM);
+  // the following cuts are limit the objects used to calculate the image zero points.
+  // These are used in bcatalog.c and apply to Mcat (nominal measurement magnitude)
+  GetConfig (config, "SIGMA_LIM",              "%lf",    0, &SIGMA_LIM);
+  GetConfig (config, "MAG_LIM",                "%lf",    0, &MAG_LIM);
+
+  DefConfig ("MAG_MIN",                        "%lf", 15.0, MAG_MIN);
+  DefConfig ("STAR_MAX_PSF_KRON",              "%lf",  0.2, STAR_MAX_PSF_KRON);
+  DefConfig ("PSF_QF_MIN",                     "%lf", 0.95, PSF_QF_MIN);
 
   DefConfig ("STAR_SCATTER",                    "%lf",  0.1, STAR_SCATTER);

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

@@ -531 +531 @@
 
   // until the analysis has converged a bit, do not use the IRLS analysis
-  // default is MaxIterations = 10
-  if (UseStandardOLS(ZPT_IMAGES)) {
+  // default is MaxIterations = 10.
+
+  // However, for stack and warp calibration, we are not iteratively solving for the
+  // zero points.  In this case, use IRLS in a single pass
+  if (!CALIBRATE_STACKS_AND_WARPS && UseStandardOLS(ZPT_IMAGES)) {
     brightStars.MaxIterations = 0;
     kronStars.MaxIterations = 0;
@@ -672 +675 @@
 
       if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) {
-        fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f  %6.3f  %6.3f  %6.3f  %6.3f : %6.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysKron, MrelKron, Mmos, Mgrid, Mflat, kronStars.alldata->yVector[Nkron]);
+        fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f %6.3f | %6.3f %6.3f | %6.3f %6.3f %6.3f %6.3f | %5.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, MsysKron, MrelKron, Mmos, Mgrid, Mflat, Moff, psfStars.alldata->dyVector[Nref]);
       }
 
@@ -702 +705 @@
     if (mark) continue;
 
-    // use liststats to find the 20-pct, median, 80-pct points
-    StatType stats;
-    liststats_setmode (&stats, "MEDIAN");
-    liststats (psfStars.alldata->yVector, NULL, NULL, Nref, &stats);
-    double altSigma = (stats.Upper80 - stats.Lower20) / 1.6;  // 20% to 80% encompasses 60% of the values, corresponds to the range (-0.85 sigma : +0.85 sigma)
-
-    // soften the individual errors with 10% of the scatter above
-    for (j = 0; j < Nref; j++) {
-      double newSigma = hypot(psfStars.alldata->dyVector[j], 0.1*altSigma);
-      psfStars.alldata->dyVector[j] = newSigma;
-    }
-
-    // soften the errors based on the scatter
+    // soften the errors based on the core scatter
     FitDataSetSoften (&psfStars, Nref);
 

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

@@ -34 +34 @@
 }  
 
-// 
+// Nsecfilt is a property of the photcode table.  we only need to retrieve it once.
+static int Nsecfilt = -1;
+
+// return the average photometry for this measurement
 float getMrel (Catalog *catalog, off_t meas, int cat, dvoMagClassType class, dvoMagSourceType source) {
 
@@ -41 +44 @@
 
   int ecode = GetPhotcodeEquivCodebyCode (photcode);
+  if (USE_REF_EQUIV) {
+    // tie this photometry to an alternative refernce
+    // (e.g., tie UDR4.i to PV3.i)
+    int refcode = GetPhotcodeEquivCodebyCode (ecode);
+    if (refcode < 0) return NAN;
+    ecode = refcode;
+  }
+
   int Nsec = GetPhotcodeNsec(ecode);
-  int Nsecfilt = GetPhotcodeNsecfilt ();
+  if (Nsec < 0) return NAN;
+
+  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
 
   int entry = Nsecfilt*ave+Nsec;
@@ -312 +325 @@
   }
 
-  int Nsecfilt = GetPhotcodeNsecfilt ();
+  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
 
   SetMrelInfo summary, results;
@@ -334 +347 @@
   int i;
 
-  int Nsecfilt = GetPhotcodeNsecfilt ();
+  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
 
   SetMrelInfo summary, results;
@@ -424 +437 @@
   ThreadInfo *threadinfo = data;
 
-  int Nsecfilt = GetPhotcodeNsecfilt ();
+  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
 
   SetMrelInfo results;
@@ -451 +464 @@
   off_t k;
 
-  int Nsecfilt = GetPhotcodeNsecfilt ();
+  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
 
   off_t m = average[0].measureOffset;
@@ -559 +572 @@
   off_t j;
 
-  int Nsecfilt = GetPhotcodeNsecfilt ();
+  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
 
   for (i = 0; i < Ncatalog; i++) {
@@ -601 +614 @@
   ALLOCATE (slist, double, Ntot);
 
-  int Nsecfilt = GetPhotcodeNsecfilt ();
+  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
 
   // eliminate bad stars using the stats for a single secfilt at a time
@@ -666 +679 @@
   // N1 = N2 = N3 = N4 = N0 = 0;
 
-  int Nsecfilt = GetPhotcodeNsecfilt ();
+  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
 
   Ntot = 0;
@@ -719 +732 @@
   StatType stats;
 
-  int Nsecfilt = GetPhotcodeNsecfilt ();
+  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
 
   Ntot = 0;
@@ -760 +773 @@
   StatType stats;
 
-  int Nsecfilt = GetPhotcodeNsecfilt ();
+  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
 
   Ntot = 0;
@@ -805 +818 @@
   ALLOCATE (Mlist, double, NBIN);
 
-  int Nsecfilt = GetPhotcodeNsecfilt ();
+  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
 
   int Ns;
@@ -840 +853 @@
   Graphdata graphdata;
 
-  int Nsecfilt = GetPhotcodeNsecfilt ();
+  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
 
   Ntotal = 0;
@@ -911 +924 @@
   return;
 
-  int Nsecfilt = GetPhotcodeNsecfilt ();
+  if (Nsecfilt < 0) { Nsecfilt = GetPhotcodeNsecfilt (); }
 
   FILE *fout = NULL;

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

@@ -237 +237 @@
   }
 
+  USE_REF_EQUIV = FALSE;
+  if ((N = get_argument (argc, argv, "-use-ref-equiv"))) {
+    remove_argument (N, &argc, argv);
+    USE_REF_EQUIV = TRUE;
+  }
   PRESERVE_PS1 = FALSE;
   if ((N = get_argument (argc, argv, "-preserve-ps1"))) {
@@ -521 +526 @@
     remove_argument (N, &argc, argv);
     USE_REFERENCE_WEIGHT = TRUE;
+  }
+
+  /* Normally, stacks and warps are treated differently from chips.  In this case, 
+     WHAT IS DIFFERENT?
+   */
+
+  CALIBRATE_WARPS_LIKE_CHIPS = FALSE;
+  if ((N = get_argument (argc, argv, "-calibrate-warps-like-chips"))) {
+    remove_argument (N, &argc, argv);
+    CALIBRATE_WARPS_LIKE_CHIPS = TRUE;
+  }
+  CALIBRATE_STACKS_LIKE_CHIPS = FALSE;
+  if ((N = get_argument (argc, argv, "-calibrate-stacks-like-chips"))) {
+    remove_argument (N, &argc, argv);
+    CALIBRATE_STACKS_LIKE_CHIPS = TRUE;
   }
 
@@ -836 +856 @@
   }
 
+  USE_REF_EQUIV = FALSE;
+  if ((N = get_argument (argc, argv, "-use-ref-equiv"))) {
+    remove_argument (N, &argc, argv);
+    USE_REF_EQUIV = TRUE;
+  }
   PRESERVE_PS1 = FALSE;
   if ((N = get_argument (argc, argv, "-preserve-ps1"))) {
@@ -978 +1003 @@
   }
 
+  CALIBRATE_WARPS_LIKE_CHIPS = FALSE;
+  if ((N = get_argument (argc, argv, "-calibrate-warps-like-chips"))) {
+    remove_argument (N, &argc, argv);
+    CALIBRATE_WARPS_LIKE_CHIPS = TRUE;
+  }
+  CALIBRATE_STACKS_LIKE_CHIPS = FALSE;
+  if ((N = get_argument (argc, argv, "-calibrate-stacks-like-chips"))) {
+    remove_argument (N, &argc, argv);
+    CALIBRATE_STACKS_LIKE_CHIPS = TRUE;
+  }
+
   STAGES = STAGE_CHIP | STAGE_WARP | STAGE_STACK;
   if ((N = get_argument (argc, argv, "-skip-chip"))) {

trunk/Ohana/src/relphot/src/bcatalog.c

@@ -1 +1 @@
 # include "relphot.h"
+
+# define BAD_COUNT(FIELD) \
+  badCount.FIELD ++; \
+  if (doTest1 && (catalog[0].measure[offset].imageID == TEST_IMAGE1)) { badTest1.FIELD ++; } \
+  if (doTest2 && (catalog[0].measure[offset].imageID == TEST_IMAGE2)) { badTest2.FIELD ++; }
+
+typedef struct {
+  int Ncode;
+  int Ntime;
+  int Ndophot;
+  int Nmag;
+  int Nsigma;
+  int Nimag;
+  int Nfew;
+  int Ngalaxy;
+  int Npsfqf;
+  int Nnan;
+  int Nbad;
+  int Npoor;
+} MyBadCountType;
 
 int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog);
@@ -9 +29 @@
   off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
   float mag;
-  int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf, Nnan, Nbad, Npoor;
+
+  // int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf, Nnan, Nbad, Npoor;
 
   int Nsecfilt = GetPhotcodeNsecfilt ();
@@ -26 +47 @@
   int NmMin = catalog[0].Nmeasure;
 
-  Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = Nnan = Nbad = Npoor = 0;
+  MyBadCountType badCount = {0,0,0, 0,0,0, 0,0,0, 0,0,0};
+  MyBadCountType badTest1 = {0,0,0, 0,0,0, 0,0,0, 0,0,0};
+  MyBadCountType badTest2 = {0,0,0, 0,0,0, 0,0,0, 0,0,0};
+
+  int doTest1 = (TEST_IMAGE1 >= 0);
+  int doTest2 = (TEST_IMAGE2 >= 0);
+
+  FILE *ftest1 = (doTest1) ? fopen ("test1.meas.dat", "a") : NULL;
+  FILE *ftest2 = (doTest2) ? fopen ("test2.meas.dat", "a") : NULL;
+
+  // Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = Nnan = Nbad = Npoor = 0;
 
   // copy the following fields to the subcatalog:
@@ -57 +88 @@
 
       offset = catalog[0].average[i].measureOffset + j;
+
+      if (doTest1) {
+        Measure *m = &catalog[0].measure[offset];
+        if (m->imageID == TEST_IMAGE1) {
+          fprintf (ftest1, "%5d %5d %5d | %5d %12d %5.3f %5.3f %5.3f %5.3f %4.3f %4.3f 0x%08x\n", m->detID, m->objID, m->catID, m->photcode, m->t, m->M, m->Map, m->Mkron, m->dM, m->psfQF, m->psfQFperf, m->photFlags);
+        }
+      }
+      if (doTest2) {
+        Measure *m = &catalog[0].measure[offset];
+        if (m->imageID == TEST_IMAGE2) {
+          fprintf (ftest2, "%5d %5d %5d | %5d %12d %5.3f %5.3f %5.3f %5.3f %4.3f %4.3f 0x%08x\n", m->detID, m->objID, m->catID, m->photcode, m->t, m->M, m->Map, m->Mkron, m->dM, m->psfQF, m->psfQFperf, m->photFlags);
+        }
+      }
 
       /* select measurements by photcode */
@@ -69 +113 @@
       }
       if (!found) {
-        Ncode ++; 
+        BAD_COUNT(Ncode);
         continue; 
       }
@@ -75 +119 @@
       /* select measurements by time */
       if (TimeSelect) {
-        if (catalog[0].measure[offset].t < TSTART) { Ntime ++; continue; }
-        if (catalog[0].measure[offset].t > TSTOP)  { Ntime ++; continue; }
+        if (catalog[0].measure[offset].t < TSTART) { BAD_COUNT(Ntime); continue; }
+        if (catalog[0].measure[offset].t > TSTOP)  { BAD_COUNT(Ntime); continue; }
       }
 
       /* select measurements by quality */
-      if (DophotSelect && ((catalog[0].measure[offset].photFlags >> 16) != DophotValue)) { Ndophot ++; continue; }
+      if (DophotSelect && ((catalog[0].measure[offset].photFlags >> 16) != DophotValue)) { BAD_COUNT(Ndophot); continue; }
 
       float Maper = USE_APER_FOR_STARGAL ? catalog[0].measure[offset].Map : catalog[0].measure[offset].Mkron;
 
       // skip garbage measurements
-      if (isnan(catalog[0].measure[offset].psfQF)     || (catalog[0].measure[offset].psfQF < 0.95))     { Npsfqf ++; continue; }
-      if (isnan(catalog[0].measure[offset].psfQFperf) || (catalog[0].measure[offset].psfQFperf < 0.95)) { Npsfqf ++; continue; }
-      if (isnan(catalog[0].measure[offset].M)) { Nnan ++; continue; }
-      if (isnan(Maper)) { Nnan ++; continue; }
+      if (isnan(catalog[0].measure[offset].psfQF)     || (catalog[0].measure[offset].psfQF     < PSF_QF_MIN)) { BAD_COUNT(Npsfqf); continue; }
+      if (isnan(catalog[0].measure[offset].psfQFperf) || (catalog[0].measure[offset].psfQFperf < PSF_QF_MIN)) { BAD_COUNT(Npsfqf); continue; }
+      if (isnan(catalog[0].measure[offset].M)) { BAD_COUNT(Nnan); continue; }
+      if (isnan(Maper)) { BAD_COUNT(Nnan); continue; }
 
       // require 0x01 in photFlags (fitted with a PSF) -- add to the test data
-      if (REQUIRE_PSFFIT && (catalog[0].measure[offset].photFlags & 0x01) == 0) { Nbad ++; continue; }
+      if (REQUIRE_PSFFIT && (catalog[0].measure[offset].photFlags & 0x01) == 0) { BAD_COUNT(Nbad); continue; }
 
       // very loose cut on PSF - APER mag (Map or Mkron) -- this lightly rejects galaxies & other oddities
       float Mkp = catalog[0].measure[offset].M - Maper;
-      if (fabs(Mkp) > 1.0) { Nbad ++; continue; }
+      if (fabs(Mkp) > 1.0) { BAD_COUNT(Nbad); continue; }
 
       if (catalog[0].measure[offset].photFlags & code->photomBadMask) { 
-        Nbad++; 
+        BAD_COUNT(Nbad); 
         continue;
       }
-      if (catalog[0].measure[offset].photFlags & code->photomPoorMask) { Npoor++; continue;}
+      if (catalog[0].measure[offset].photFlags & code->photomPoorMask) { BAD_COUNT(Npoor); continue;}
 
       // weak test for galaxies
-      if (Mkp > 0.5) {
+      if (Mkp > STAR_MAX_PSF_KRON) {
         nEXT ++;
       } else {
@@ -112 +156 @@
       /* select measurements by mag limit */
       mag = PhotCat (&catalog[0].measure[offset], MAG_CLASS_PSF);
-      if (mag > MAG_LIM) { Nmag ++; continue; } 
+      if (mag > MAG_LIM) { BAD_COUNT(Nmag); continue; } 
+      if (mag < MAG_MIN) { BAD_COUNT(Nmag); continue; } 
+      // note that these two cuts are applied to all photcodes. these are also applied to
+      // the nominal apparent mags : if the expected sensitivity is very different between
+      // different filters, or if exposure times have a large range, these cuts may not be
+      // effective (or may exclude too much).
 
       /* select measurements by measurement error */
-      if ((SIGMA_LIM > 0) && (catalog[0].measure[offset].dM > SIGMA_LIM)) { Nsigma ++; continue; }
+      if ((SIGMA_LIM > 0) && (catalog[0].measure[offset].dM > SIGMA_LIM)) { BAD_COUNT(Nsigma); continue; }
 
       /* select measurements by mag limit */
       if (ImagSelect) {
         mag = PhotInst (&catalog[0].measure[offset], MAG_CLASS_PSF);
-        if (mag < ImagMin) { Nimag ++; continue; }
-        if (mag > ImagMax) { Nimag ++; continue; }
-      }
+        if (mag < ImagMin) { BAD_COUNT(Nimag); continue; }
+        if (mag > ImagMax) { BAD_COUNT(Nimag); continue; }
+      }
+      // note that instrumental mag cuts are reasonable for individual exposures using the same
+      // detector, but are probably not appropriate for stacks with a large range of the number
+      // of inputs.  
 
       // count this measurement as valid for this secfilt entry
@@ -147 +199 @@
     if (nEXT >= nPSF) {
       Nmeasure -= Nm;
-      Ngalaxy ++;
+      badCount.Ngalaxy++;
       continue; 
     }
@@ -157 +209 @@
     if (Nm <= STAR_TOOFEW) { /* enough measurements total? */
       Nmeasure -= Nm;
-      Nfew ++;
+      badCount.Nfew++;
       continue; 
     }
@@ -167 +219 @@
     if (!anySecfiltGood) {
       Nmeasure -= Nm;
-      Nfew ++;
+      badCount.Nfew++;
       continue; 
     }
@@ -204 +256 @@
              subcatalog[0].Naverage,  subcatalog[0].Nmeasure,  catalog[0].Naverage, catalog[0].Nmeasure, catalog[0].filename, NmMin, NmMax);
     fprintf (stderr, "rejections: %d stars have too few measures:\n   %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n", 
-             Nfew, Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Npsfqf, Nnan, Ngalaxy, Nbad, Npoor);
+             badCount.Nfew, badCount.Ncode, badCount.Ntime, badCount.Ndophot, badCount.Nmag, badCount.Nsigma, badCount.Nimag, badCount.Npsfqf, badCount.Nnan, badCount.Ngalaxy, badCount.Nbad, badCount.Npoor);
+  }
+
+  if (doTest1) {
+    fprintf (stderr, "TEST1 rejections: %d stars have too few measures: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n", 
+             badTest1.Nfew, badTest1.Ncode, badTest1.Ntime, badTest1.Ndophot, badTest1.Nmag, badTest1.Nsigma, badTest1.Nimag, badTest1.Npsfqf, badTest1.Nnan, badTest1.Ngalaxy, badTest1.Nbad, badTest1.Npoor);
+  }
+  if (doTest2) {
+    fprintf (stderr, "TEST2 rejections: %d stars have too few measures: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n", 
+             badTest2.Nfew, badTest2.Ncode, badTest2.Ntime, badTest2.Ndophot, badTest2.Nmag, badTest2.Nsigma, badTest2.Nimag, badTest2.Npsfqf, badTest2.Nnan, badTest2.Ngalaxy, badTest2.Nbad, badTest2.Npoor);
   }
 
@@ -212 +273 @@
   }
 
-  free (Nvalid)
+  free (Nvalid);
+
+  if (ftest1) { fclose (ftest1); }
+  if (ftest2) { fclose (ftest2); }
 
   return (TRUE);

trunk/Ohana/src/relphot/src/extra.c

@@ -70 +70 @@
 }
 
-// for now (20140710) I need to identify gpc1 stacks explicitly.  generalize in the future
+// for now (20140710) I need to identify gpc1 & gpc2 stacks explicitly.  generalize in the future
 int isGPC1stack (int photcode) {
 
-  if ((photcode == 11000) || (photcode == 14100)) return TRUE; // g-band
-  if ((photcode == 11100) || (photcode == 14200)) return TRUE; // r-band
-  if ((photcode == 11200) || (photcode == 14300)) return TRUE; // i-band
-  if ((photcode == 11300) || (photcode == 14400)) return TRUE; // z-band
-  if ((photcode == 11400) || (photcode == 14500)) return TRUE; // y-band
-  if ((photcode == 11500) || (photcode == 14600)) return TRUE; // w-band
+  // 11000 : GPC1.g.SkyChip
+  // 11001 : VARIOUS.g.SkyChip (output of xcamera stacks)
+  // 14100 : SIMTEST.g.SkyChip
+  // 31000 : GPC2.g.SkyChip
+  
+  if (CALIBRATE_STACKS_LIKE_CHIPS) return FALSE;
+
+  if ((photcode == 11000) || (photcode == 11001) || (photcode == 14100) || (photcode == 31000)) return TRUE; // g-band
+  if ((photcode == 11100) || (photcode == 11101) || (photcode == 14200) || (photcode == 31100)) return TRUE; // r-band
+  if ((photcode == 11200) || (photcode == 11201) || (photcode == 14300) || (photcode == 31200)) return TRUE; // i-band
+  if ((photcode == 11300) || (photcode == 11301) || (photcode == 14400) || (photcode == 31300)) return TRUE; // z-band
+  if ((photcode == 11400) || (photcode == 11401) || (photcode == 14500) || (photcode == 31400)) return TRUE; // y-band
+  if ((photcode == 11500) || (photcode == 11501) || (photcode == 14600) || (photcode == 31500)) return TRUE; // w-band
 
   return FALSE;
@@ -85 +92 @@
 // for now (20140710) I need to identify gpc1 stacks explicitly.  generalize in the future
 int isGPC1warp (int photcode) {
+
+  if (CALIBRATE_WARPS_LIKE_CHIPS) return FALSE;
+
+  // 11000 : GPC1.g.SkyChip
+  // 14100 : SIMTEST.g.SkyChip
+
+  // 12000 : GPC1.g.ForcedWarp
+  // 15100 : SIMTEST.g.ForcedWarp
+  // 32000 : GPC2.g.ForcedWarp
 
   // diff warps get stack-like photcodes (kind of lame)
@@ -97 +113 @@
   }
 
-  if ((photcode == 12000) || (photcode == 15100)) return TRUE; // g-band
-  if ((photcode == 12100) || (photcode == 15200)) return TRUE; // r-band
-  if ((photcode == 12200) || (photcode == 15300)) return TRUE; // i-band
-  if ((photcode == 12300) || (photcode == 15400)) return TRUE; // z-band
-  if ((photcode == 12400) || (photcode == 15500)) return TRUE; // y-band
-  if ((photcode == 12500) || (photcode == 15600)) return TRUE; // w-band
+  if ((photcode == 12000) || (photcode == 15100) || (photcode == 32000)) return TRUE; // g-band
+  if ((photcode == 12100) || (photcode == 15200) || (photcode == 32100)) return TRUE; // r-band
+  if ((photcode == 12200) || (photcode == 15300) || (photcode == 32200)) return TRUE; // i-band
+  if ((photcode == 12300) || (photcode == 15400) || (photcode == 32300)) return TRUE; // z-band
+  if ((photcode == 12400) || (photcode == 15500) || (photcode == 32400)) return TRUE; // y-band
+  if ((photcode == 12500) || (photcode == 15600) || (photcode == 32500)) return TRUE; // w-band
 
   return FALSE;

trunk/Ohana/src/relphot/src/fit1d_irls.c

@@ -374 +374 @@
   StatType stats;
   liststats_setmode (&stats, "MEDIAN");
-  liststats (dataset->alldata->yVector, NULL, NULL, Nvalues, &stats);
+  liststats (dataset->alldata->yVector, dataset->alldata->dyVector, NULL, Nvalues, &stats);
+  // liststats sorts yVector and dyVector together
 
   double altSigma = (stats.Upper80 - stats.Lower20) / 1.6;  // 20% to 80% encompasses 60% of the values, corresponds to the range (-0.85 sigma : +0.85 sigma)
@@ -380 +381 @@
   // soften the individual errors with 10% of the scatter above
   for (int j = 0; j < Nvalues; j++) {
-    double newSigma = hypot(dataset->alldata->dyVector[j], 0.1*altSigma);
+    double newSigma = hypot(dataset->alldata->dyVector[j], 0.5*altSigma);
     dataset->alldata->dyVector[j] = newSigma;
   }

trunk/Ohana/src/relphot/src/initialize.c

@@ -50 +50 @@
       fprintf (stderr, "DophotSelect: FALSE\n");
     }
-    fprintf (stderr, "PSF_QF limit: 0.85 (hardwired)\n");
+    fprintf (stderr, "PSF_QF limit: %f\n", PSF_QF_MIN);
 
     // fprintf (stderr, "Photom Bad Mask: 0x%08x, Photom Poor Mask: 0x%08x\n");
 
-    fprintf (stderr, "MAG_LIM: %f, SIGMA_LIM: %f\n", MAG_LIM, SIGMA_LIM);
+    fprintf (stderr, "MAG_MIN: %f, MAG_LIM: %f, SIGMA_LIM: %f\n", MAG_MIN, MAG_LIM, SIGMA_LIM);
     fprintf (stderr, "INST_MAG_MIN: %f, INST_MAG_MAX: %f\n", ImagMin, ImagMax);
 
+    fprintf (stderr, "STAR_MAX_PSF_KRON: %f\n", STAR_MAX_PSF_KRON);
     fprintf (stderr, "STAR_TOOFEW: %d\n", STAR_TOOFEW);
 

trunk/Ohana/src/relphot/src/load_catalogs.c

@@ -151 +151 @@
     // TimeSelect -time
     // DophotSelect
-    // (note that psfQF is applied rigidly at 0.85, as is the galaxy test)
+    // PSF_QF_MIN
+    // STAR_MAX_PSF_KRON
+    // MAG_MIN
     // MAG_LIM
     // SIGMA_LIM
@@ -158 +160 @@
 
     char *command = NULL;
-    strextend (&command, "relphot_client %s -load %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -D CAMERA %s -D MAG_LIM %f -D SIGMA_LIM %f", 
-              PhotcodeList, table->hosts[i].results, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, CAMERA, MAG_LIM, SIGMA_LIM);
+    strextend (&command, "relphot_client %s -load %s -hostID %d -D CATDIR %s -hostdir %s", PhotcodeList, table->hosts[i].results, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname);
+    strextend (&command, "-region %f %f %f %f -D CAMERA %s", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, CAMERA);
+    strextend (&command, "-D MAG_MIN %f -D MAG_LIM %f -D SIGMA_LIM %f", MAG_MIN, MAG_LIM, SIGMA_LIM);
+    strextend (&command, "-D STAR_MAX_PSF_KRON %f -D PSF_QF_MIN %f", STAR_MAX_PSF_KRON, PSF_QF_MIN);
 
     if (VERBOSE)             { strextend (&command, "-v"); }

trunk/Ohana/src/relphot/src/reload_catalogs.c

@@ -251 +251 @@
     if (!(STAGES & STAGE_WARP))  { strextend (&command, "-skip-warp"); }
     if (!(STAGES & STAGE_STACK)) { strextend (&command, "-skip-stack"); }
+    if (CALIBRATE_WARPS_LIKE_CHIPS) { strextend (&command, "-calibrate-warps-like-chips"); }
+    if (CALIBRATE_STACKS_LIKE_CHIPS) { strextend (&command, "-calibrate-stacks-like-chips"); }
 
     // deprecate

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

@@ -33 +33 @@
   
   switch (mode) {
-    case UPDATE_IMAGES:
+    case UPDATE_IMAGES: // -images
       // calculate zero points for images (may group by exposure [mosaic], night [tgroups]; may measure flat-correction)
-      // IF CALLED WITH NLOOP == 0, DOES NOT LOAD bcatalog, just loads image table and generates ImageSubset.dat
+      // IF CALLED WITH NLOOP == 0 (and NOT -only-stacks-and-warps), this step DOES NOT LOAD bcatalog.  Instead it just loads
+      // the image table and applies the Mcal values to the measurements.  In parallel context, it generates ImageSubset.dat.
       // If called with -update, calls reload_catalogs() just like apply_offsets, UNLESS -only-stacks-and-warps is selected
       relphot_images (skylist);