Index: /trunk/Ohana/src/relphot/include/relphot.h
===================================================================
--- /trunk/Ohana/src/relphot/include/relphot.h	(revision 39463)
+++ /trunk/Ohana/src/relphot/include/relphot.h	(revision 39464)
@@ -132,4 +132,8 @@
   double *psfqfperf_list;     // psfqfperf for all filters
   double *stargal_list;	      // stargal for all filters
+  uint32_t *photflag_list;      // photflags for all filters
+
+  int Nstargal;
+  int Nphotflags;
 
   int 	*havePS1;	   // this secfilt has synthetic mags
Index: /trunk/Ohana/src/relphot/src/StarOps.c
===================================================================
--- /trunk/Ohana/src/relphot/src/StarOps.c	(revision 39463)
+++ /trunk/Ohana/src/relphot/src/StarOps.c	(revision 39464)
@@ -117,4 +117,7 @@
     results->haveUbercal[i] = 0;
 
+    results->Nstargal   = 0;
+    results->Nphotflags = 0;
+
     results->tessID[i]    = 0;
     results->projID[i]    = 0;
@@ -143,4 +146,5 @@
     ALLOCATE (results->psfqfperf_list, double, Nmax);
     ALLOCATE (results->stargal_list,   double, Nmax);
+    ALLOCATE (results->photflag_list,  uint32_t, Nmax);
 
     ALLOCATE (results->havePS1,   int, results->Nsecfilt);
@@ -183,4 +187,5 @@
     results->psfqfperf_list = NULL;
     results->stargal_list   = NULL;
+    results->photflag_list  = NULL;
 
     results->havePS1 = NULL;
@@ -214,4 +219,5 @@
   FREE (results->psfqfperf_list);
   FREE (results->stargal_list);
+  FREE (results->photflag_list);
 
   FREE (results->havePS1);
@@ -423,544 +429,4 @@
     continue; }
 
-
-// setMrel_catalog is used in 3 different contexts:
-
-// * during the main relphot iterations:
-// ** operations only apply to measureT / averageT
-// ** we are applying the analysis to the bright subset catalog
-// ** we skip any stars found to be bad (STAR_BAD)
-
-// * during the final pass
-// ** operations are applied to all objects
-// ** we have special tests for PS1 extended data, 2MASS extended & good data, & synthetic photometry
-
-// * during the stand-alone '-averages' mode (relphot_objects)
-// ** no image data is loaded
-// ** we only use the provided calibration (measure.Mcal)
-// ** we skip outlier rejection of measurements
-
-// set the Mrel values for the specified catalog
-// NOTE: here 'catalog' is a pointer to a specific catalog, not the root of the array
-int setMrel_catalog (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
-  OHANA_UNUSED_PARAM(catalog);
-  OHANA_UNUSED_PARAM(Nc);
-  OHANA_UNUSED_PARAM(pass);
-  OHANA_UNUSED_PARAM(flatcorr);
-  OHANA_UNUSED_PARAM(results);
-  OHANA_UNUSED_PARAM(Nsecfilt);
-
-  fprintf (stderr, "??? should you still be using the old version of setMrel_catalog??\n");
-  exit (3);
-
-# if (0)
-
-  off_t j, k, m, ID;
-  int N;
-  float Msys, Mcal= 0, Mmos = 0, Mgrid = 0;
-
-  StatType stats, apstats, kronstats;
-  liststats_setmode (&psfstats, STATMODE);
-  liststats_setmode (&apstats, STATMODE);
-  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;
-
-  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
-
-  int isSetMrelFinal = (pass >= 0);
-
-  char *primaryCell = NULL;
-  if (isSetMrelFinal) {
-    ALLOCATE (primaryCell, char, DVO_MAX_PATH);
-  }
-
-  for (j = 0; j < catalog[Nc].Naverage; j++) {
-    // option for a test print
-    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);
-    }
-
-    if (isSetMrelFinal) {
-      // set the name of the primary skycell (this is used in a strcmp to match the skycells in stack detections)
-      // XXX : this whole function is deprecated
-      // BoundaryTreePrimaryCell(primaryCell, catalog[Nc].average[j].R, catalog[Nc].average[j].D);
-    }
-
-    int GoodPS1 = FALSE;
-    int Good2MASS = FALSE;
-    int Galaxy2MASS = FALSE;
-
-    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;
-    for (Ns = 0; Ns < Nphotcodes; Ns++) {
-
-      int thisCode = photcodes[Ns][0].code;
-      int Nsec = GetPhotcodeNsec(thisCode);
-
-      /* calculate the average mag in this SEC photcode for a single star */
-
-      /* star/photcodes already calibrated */
-      if ( isSetMrelFinal && catalog[Nc].found_t[Nsecfilt*j+Nsec]) continue;  
-      
-      // skip bad stars
-      if (!isSetMrelFinal && (catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD)) continue;
-
-      int Ncode = 0;
-      int Next = 0;
-      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
-      // closest for the subset of primary projection cells
-
-      float stackCenterOffsetMin = 1e9;
-      off_t stackCenterMeasureMin = -1;
-
-      float stackPrimaryOffsetMin = 1e9;
-      off_t stackPrimaryMeasureMin = -1;
-
-      int forceSynth = FALSE;
-      int forceSynthEntry = -1;
-
-      int haveUbercal = FALSE;
-
-      int minUbercalDist = 1000;
-    
-      int isBad = FALSE;
-      int isSuspect = FALSE;
-
-      N = 0;
-      m = catalog[Nc].averageT[j].measureOffset;
-      for (k = 0; k < catalog[Nc].averageT[j].Nmeasure; k++, m++) {
-
-	// skip measurements that do not match the current photcode
-	PhotCode *code = GetPhotcodebyCode (catalog[Nc].measureT[m].photcode);
-	if (!code) continue;
-	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); 
-
-	if (getImageEntry (m, Nc) < 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;
-	  Mcal = catalog[Nc].measureT[m].Mcal; // check that this is zero for loaded REF value
-	} else {
-	  Mcal  = getMcal  (m, Nc, flatcorr, catalog);
-	  if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
-	  Mmos  = getMmos  (m, Nc);
-	  if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
-	  Mgrid = getMgrid (m, Nc);
-	  if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
-	}
-
-	// skip some absurd values NAN, < 0.0, > 30.0
-	Msys = PhotSysTiny (&catalog[Nc].measureT[m], &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt], MAG_CLASS_PSF);
-	if (isnan(Msys)) SKIP_THIS_MEAS(Nsys);
-	if (Msys <  0.0) SKIP_THIS_MEAS(Nsys);
-	if (Msys > 30.0) SKIP_THIS_MEAS(Nsys);
-
-	int myUbercalDist = getUbercalDist(m,Nc);
-	minUbercalDist = MIN(minUbercalDist, myUbercalDist);
-
-	if (isSetMrelFinal) {
-	  // 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)) {
-	      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 ++;
-	      } else {
-		NpsfPS1 ++;
-	      }
-	    }
-	    havePS1 = TRUE;
-	  }
-
-	  // gpc1 stack data
-	  if ((catalog[Nc].measure[m].photcode >= 11000) && (catalog[Nc].measure[m].photcode <= 11400)) {
-	    haveStack = TRUE;
-	    haveStackObject = 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;
-	    }
-
-	    // 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
-	  }
-
-	  // count extended detections for 2MASS (XXX NOTE hardwired photcodes 2011, 2012, 2013)
-	  if ((catalog[Nc].measure[m].photcode >= 2011) && (catalog[Nc].measure[m].photcode <= 2013)) {
-	    if (catalog[Nc].measure[m].photFlags & 0x00c00000) {
-	      Next ++;
-	      Galaxy2MASS = TRUE;
-	    }
-	    if (pass == 0) {
-	      if (catalog[Nc].measure[m].photFlags & 0x00000007) {
-		Good2MASS = TRUE;
-	      } else {
-		// detections without one of these bits should only be used in PASS_1
-		SKIP_THIS_MEAS(Nbad);
-	      }
-	    }
-	  }
-
-	  // Blindly accepth the SYNTH mags if we are above saturation, otherwise, 
-	  // ignore SYNTH photcodes until PASS == 4 (where we also accept saturated stars)
-	  if ((catalog[Nc].measure[m].photcode >= 3001) && (catalog[Nc].measure[m].photcode <= 3005)) {
-	    // something of a hack: force object to use synth values if synth mags >>
-	    // saturation (3pi instrumental mags < -15)
-	    float MaxMagForceSynth = NAN;
-	    switch (catalog[Nc].measure[m].photcode) {
-	      case 3001:
-		MaxMagForceSynth = 13.64;
-		break;
-	      case 3002:
-		MaxMagForceSynth = 13.76;
-		break;
-	      case 3003:
-		MaxMagForceSynth = 13.74;
-		break;
-	      case 3004:
-		MaxMagForceSynth = 12.94;
-		break;
-	      case 3005:
-		MaxMagForceSynth = 12.01;
-		break;
-	    }
-	    if (catalog[Nc].measureT[m].M < MaxMagForceSynth) {
-	      forceSynth = TRUE;
-	      forceSynthEntry = N;
-	    } else {
-	      if (pass < 4) {
-		SKIP_THIS_MEAS(Nbad);
-	      }
-	      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++;
-	}
-
-	// dlist gives the error per measurement, wlist gives the weight
-	// we can modify the error and weight in a few ways:
-	// 1) MIN_ERROR guarantees a floor
-	// 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
-	// 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
-	// 4) some reference photcode of some kind can be specified as fixed and have a high weight
-	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...)
-	if (catalog[Nc].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
-	  haveUbercal = TRUE;
-	  wlist[N] = 10.0;
-	}
-
-	// tie down reference photometry if the -refcode (code) option is selected
-	// eg, -refcode g_SDSS
-	// this probably makes no sense in the context of multifilter analysis
-	if (refPhotcode) {
-	  if (code->code == refPhotcode->code) {
-	    wlist[N] = 100.0;
-	  }
-	}
-	N++;
-      }
-
-      int Nminmeas = isSetMrelFinal ? 1 : STAR_TOOFEW + 1;
-
-      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 ++;
-	}
-      }
-
-      // when performing the grid analysis, STAR_TOOFEW should be set to 1;
-      if (N < Nminmeas) { /* too few measurements */
-	// fprintf (f, "%10.6f %10.6f %d %d %d\n", catalog[Nc].averageT[j].R, catalog[Nc].averageT[j].D, catalog[Nc].measureT[catalog[Nc].averageT[j].measureOffset].imageID, N, STAR_TOOFEW); 
-	catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_OBJ_FEW;
-	if (Ncode == 0) {
-	  results->Ncode ++;
-	} else {
-	  results->Nfew ++;
-	}
-	continue;
-      } else {
-	catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_OBJ_FEW;
-      }	
-
-      if (forceSynth) {
-	// use the single SYNTH value instead of the other mags here
-	myAssert ((forceSynthEntry < N) && (forceSynthEntry >= 0), "programming error");
-	list[0]  = list[forceSynthEntry];
-	dlist[0] = dlist[forceSynthEntry];
-	wlist[0] = wlist[forceSynthEntry];
-	N = 1;
-      }
-      liststats (list, dlist, wlist, N, &stats);
-
-      catalog[Nc].secfilt[Nsecfilt*j+Nsec].M      = stats.mean;
-      catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM     = stats.error;
-      catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mchisq = (stats.Nmeas > 1) ? stats.chisq : NAN;
-
-      // when running -averages, we have no information about the images, so we cannot set this
-      if (minUbercalDist > -1) {
-	catalog[Nc].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist;
-      }
-
-      if (isSetMrelFinal) {
-	catalog[Nc].found_t[Nsecfilt*j+Nsec] = TRUE;
-
-	catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mstdev = 1000.0*stats.sigma; // Mstdev is in millimags (not enough space for more precision)
-	catalog[Nc].secfilt[Nsecfilt*j+Nsec].Ncode = Ncode;
-	catalog[Nc].secfilt[Nsecfilt*j+Nsec].Nused = stats.Nmeas;
-
-	catalog[Nc].secfilt[Nsecfilt*j+Nsec].M_80 = 1000 * stats.Upper80;
-	catalog[Nc].secfilt[Nsecfilt*j+Nsec].M_20 = 1000 * stats.Lower20;
-
-	// NOTE : use the modified weight for apmags as well as psf mags
-	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;
-
-	  // 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 = 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
-	  // 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].stackDetectID = ID;
-
-	  catalog[Nc].measure[m].dbFlags |= ID_MEAS_STACK_PHOT_SRC;
-	  if (stackPrimaryMeasureMin >= 0) {
-	    catalog[Nc].measure[m].dbFlags |= ID_MEAS_STACK_PRIMARY;
-	  }
-	}
-
-	// NOTE: for 2MASS measurements, Next should be 1, as should N
-	if ((Next > 0) && (Next > 0.5*N)) {
-	  catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_OBJ_EXT;
-	}
-
-	switch (pass) {
-	  case 0:
-	    catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_0;
-	    if (havePS1) GoodPS1 = TRUE;
-	    break;
-	  case 1:
-	    catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_1;
-	    if (havePS1) GoodPS1 = TRUE;
-	    break;
-	  case 2:
-	    catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_2;
-	    if (havePS1) GoodPS1 = TRUE;
-	    break;
-	  case 3:
-	    catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_3;
-	    break;
-	  case 4:
-	    catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_4;
-	    break;
-	}
-	if (haveSynth) {
-	  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;
-	}	
-      }
-    }
-
-    if (isSetMrelFinal) {
-      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)) {
-	catalog[Nc].average[j].flags |= ID_OBJ_EXT;
-      }
-      if (GoodPS1) {
-	catalog[Nc].average[j].flags |= ID_OBJ_GOOD;
-      }
-      if (Galaxy2MASS) {
-	catalog[Nc].average[j].flags |= ID_OBJ_EXT_ALT;
-      }
-      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;
-      }
-    }
-  }
-  if (primaryCell) free (primaryCell);
-  return (TRUE);
-# endif
-}
 
 int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {
Index: /trunk/Ohana/src/relphot/src/relphot_objects.c
===================================================================
--- /trunk/Ohana/src/relphot/src/relphot_objects.c	(revision 39463)
+++ /trunk/Ohana/src/relphot/src/relphot_objects.c	(revision 39464)
@@ -74,8 +74,11 @@
       
       for (j = 0; j < catalog.Naverage; j++) {
-	catalog.average[j].flags &= ~photomBits; // reset all except astrometry bits
+	catalog.average[j].flags    &= ~photomBits; // reset all except astrometry bits
 	catalog.average[j].psfQF     = NAN;     // reset (will be re-calculated here)
 	catalog.average[j].psfQFperf = NAN;	// reset (will be re-calculated here)
 	catalog.average[j].stargal   = NAN;	// reset (will be re-calculated here)
+	catalog.average[j].photFlagsUpper = 0;	// reset (will be re-calculated here)
+	catalog.average[j].photFlagsLower = 0;	// reset (will be re-calculated here)
+	catalog.average[j].NwarpOK        = 0;	// reset (will be re-calculated here)
 	for (k = 0; k < Nsecfilt; k++) {
 	  dvo_secfilt_init (&catalog.secfilt[j*Nsecfilt + k], SECFILT_RESET_ALL);
Index: /trunk/Ohana/src/relphot/src/setMrelCatalog.c
===================================================================
--- /trunk/Ohana/src/relphot/src/setMrelCatalog.c	(revision 39463)
+++ /trunk/Ohana/src/relphot/src/setMrelCatalog.c	(revision 39464)
@@ -9,4 +9,5 @@
 
 int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags flags);
+void GetPhotFlagStats (uint32_t *photFlagUpper, uint32_t *photFlagLower, uint32_t *photflag_list, int Nphotflag);
 
 # define UBERCAL_WEIGHT 100.0
@@ -183,5 +184,11 @@
 	if (isfinite(measure[k].extNsigma) && isfinite(measure[k].psfQF) && (measure[k].psfQF > 0.85)) {
 	  stargalmax = MAX (stargalmax, measure[k].extNsigma);
+	  results->stargal_list[results->Nstargal] = measure[k].extNsigma;
+	  results->Nstargal ++;
 	}
+
+	// for all GPC1 measurements, track the photflags
+	results->photflag_list[results->Nphotflags] = measure[k].photFlags;
+	results->Nphotflags ++;
       }
 
@@ -454,4 +461,19 @@
       average[0].flags |= ID_OBJ_GOOD_ALT;
     }
+
+    GetPhotFlagStats (&average[0].photFlagsUpper, &average[0].photFlagsLower, results->photflag_list, results->Nphotflags);
+
+    // calculate the median of the stargal (ExtNsigma) values
+    dsort(results->stargal_list, results->Nstargal);
+    double stargal_median;
+    if (results->Nstargal % 2) {
+      int Ncenter = results->Nstargal / 2;
+      stargal_median = results->stargal_list[Ncenter];
+    } else {
+      int Ncenter = results->Nstargal / 2 - 1;
+      stargal_median = 0.5*(results->stargal_list[Ncenter] + results->stargal_list[Ncenter + 1]);
+    }
+    average[0].stargal   = stargal_median;
+
     average[0].psfQF     = psfQfMax;
     average[0].psfQFperf = psfQfPerfMax;
@@ -756,4 +778,6 @@
   }
 
+  average->NwarpOK = 0; // count the number of warps with psfQF > 0.0
+
   // assign measurements to the photcode lists 
   for (k = 0; k < average->Nmeasure; k++) {
@@ -774,5 +798,6 @@
 
     results->Nmeas[Nsec] ++;
-    if (measure[k].psfQFperf > 0.85) results->NmeasGood[Nsec] ++;
+    if (isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > 0.85)) results->NmeasGood[Nsec] ++;
+    if (isfinite(measure[k].psfQF)     && (measure[k].psfQF     > 0.00)) average->NwarpOK ++;
 
     // use primary skycell for DIFF, stack skycell for WARP
@@ -960,2 +985,37 @@
 }
 
+void GetPhotFlagStats (uint32_t *photFlagUpper, uint32_t *photFlagLower, uint32_t *photflag_list, int Nphotflag) {
+
+  int i, bit;
+
+  uint32_t outputUpper = 0;
+  uint32_t outputLower = 0;
+
+  uint32_t flag = 0x1; // start with the flag on the 0 bit
+  for (bit = 0; bit < 32; bit++) {
+
+    int Nraised = 0;
+    for (i = 0; i < Nphotflag; i++) {
+      if (photflag_list[i] & flag) Nraised ++;
+    }
+
+    float fraction = Nraised / (float) Nphotflag;
+
+    if (fraction < 0.25) {
+      // do nothing (no bits raised)
+    }
+    if ((fraction >= 0.25) && (fraction < 0.50)) {
+      outputLower |= flag;
+    }
+    if ((fraction >= 0.50) && (fraction < 0.75)) {
+      outputUpper |= flag;
+    }
+    if (fraction >= 0.75) {
+      outputLower |= flag;
+      outputUpper |= flag;
+    }
+    flag = (flag << 1);
+  }
+  *photFlagUpper = outputUpper;
+  *photFlagLower = outputLower;
+}
Index: /trunk/Ohana/src/relphot/src/setMrelFinal.c
===================================================================
--- /trunk/Ohana/src/relphot/src/setMrelFinal.c	(revision 39463)
+++ /trunk/Ohana/src/relphot/src/setMrelFinal.c	(revision 39464)
@@ -94,4 +94,7 @@
     catalog[0].average[i].psfQFperf = NAN;	// force recalculation below
     catalog[0].average[i].stargal   = NAN;	// force recalculation below
+    catalog[0].average[i].photFlagsUpper = 0;	// reset (will be re-calculated here)
+    catalog[0].average[i].photFlagsLower = 0;	// reset (will be re-calculated here)
+    catalog[0].average[i].NwarpOK        = 0;	// reset (will be re-calculated here)
   }