Index: /branches/eam_branches/ipp-20150419/Ohana/src/libdvo/include/dvo.h =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/libdvo/include/dvo.h (revision 38292) +++ /branches/eam_branches/ipp-20150419/Ohana/src/libdvo/include/dvo.h (revision 38293) @@ -117,4 +117,7 @@ ID_MEAS_ICRF_QSO = 0x00040000, // this measurement is an ICRF reference position ID_MEAS_IMAGE_EPOCH = 0x00080000, // this measurement is registered to the image epoch (not tied to ref catalog epoch) + ID_MEAS_PHOTOM_PSF = 0x00100000, // this measurement is used for the mean psf mag + ID_MEAS_PHOTOM_APER = 0x00200000, // this measurement is used for the mean ap mag + ID_MEAS_PHOTOM_KRON = 0x00400000, // this measurement is used for the mean kron mag } DVOMeasureFlags; @@ -786,4 +789,6 @@ off_t *found_t; off_t *foundWarp_t; + + char *measureRank; // off_t *image_t; // off_t *mosaic_t; Index: /branches/eam_branches/ipp-20150419/Ohana/src/relphot/include/relphot.h =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/relphot/include/relphot.h (revision 38292) +++ /branches/eam_branches/ipp-20150419/Ohana/src/relphot/include/relphot.h (revision 38293) @@ -99,4 +99,7 @@ double *errlist; // mag errors for a star double *wgtlist; // weights to use for mean mags + int *ranking; // weights to use for mean mags + int *measSeq; // weights to use for mean mags + int Nlist; } StatDataSet; @@ -114,4 +117,6 @@ // the final output steps. + int Nsecfilt; + StatDataSet *psfData; // one is allocated for each primary (average) photcode StatDataSet *aperData; @@ -121,4 +126,22 @@ double *psfqfperf_list; // psfqfperf for all filters double *stargal_list; // stargal for all filters + + int *havePS1; // this secfilt has synthetic mags + int *haveSYN; // this secfilt has synthetic mags + int *measSYN; // this measurement is the synthetic mag for this secfilt + int *needSYN; // this secfilt mag should use synthetic mags + float *minSYN; // minimum synthetic mag below which synthetic should be forced + + float *psfQfMax; // max psfQf value for this secfilt + float *psfQfPerfMax; // max psfQfperf value for this secfilt + + int *Next; // count of PS1 exposure (chip) measurements for this secfilt + int *NexpPS1; // count of PS1 exposure (chip) measurements for this secfilt + int *haveUbercal; // does this secfilt have any ubercal data? + + float *minUbercalDist; + + int NextPS1; // count of PS1 exposure (chip) measurements for this secfilt + int NpsfPS1; // count of PS1 exposure (chip) measurements for this secfilt StatType psfstats; @@ -473,5 +496,5 @@ int print_measure_set_alt (Average *average, SecFilt *secfilt, Measure *measure); int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt); -int setMrelAverageExposure (off_t meas, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found); +int setMrelAverageExposure (off_t measureOffset, int cat, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT); int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t Nmeasure); int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found, off_t Nmeasure); @@ -480,4 +503,10 @@ void SetMrelInfoInit (SetMrelInfo *results, int allocLists); +void SetMrelInfoFree (SetMrelInfo *results); +void SetMrelInfoReset (SetMrelInfo *results); +void SetMrelInfoResetObject (SetMrelInfo *results); + +void StatDataSetFree (StatDataSet *dataset, int Nsecfilt); +StatDataSet *StatDataSetAlloc (int Nsecfilt, int Nmax); int init_synthetic_mags (); @@ -532,5 +561,5 @@ int isGPC1synth (int photcode); int whichGPC1filter (int photcode); - +int is2MASS (int photcode); SynthZeroPoints *SynthZeroPointsLoad (char *filename); Index: /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/ImageOps.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/ImageOps.c (revision 38292) +++ /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/ImageOps.c (revision 38293) @@ -397,4 +397,15 @@ i = MeasureToImage[cat][meas]; return (i); +} + +int getImageFlags (off_t meas, int cat) { + + off_t i; + + if (!MeasureToImage) return 0; + + i = MeasureToImage[cat][meas]; + + return (image[i].flags); } Index: /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/MosaicOps.c (revision 38292) +++ /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/MosaicOps.c (revision 38293) @@ -882,4 +882,15 @@ } +int getMosaicFlags (off_t meas, int cat) { + + if (!MOSAIC_ZEROPT) return (0); + + // unassigned measurements belong to simple chips + int i = MeasureToMosaic[cat][meas]; + if (i == -1) return (0); + + return (mosaic[i].flags); +} + typedef struct { int Nfew; Index: /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/StarOps.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/StarOps.c (revision 38292) +++ /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/StarOps.c (revision 38293) @@ -64,9 +64,11 @@ ALLOCATE (dataset[i].wgtlist, double, Nmax); ALLOCATE (dataset[i].errlist, double, Nmax); + ALLOCATE (dataset[i].ranking, int, Nmax); + ALLOCATE (dataset[i].measSeq, int, Nmax); } return dataset; } -void StatDataSetFree (StatDataSet *, int Nsecfilt) { +void StatDataSetFree (StatDataSet *dataset, int Nsecfilt) { int i; @@ -76,4 +78,6 @@ FREE (dataset[i].wgtlist); FREE (dataset[i].errlist); + FREE (dataset[i].ranking); + FREE (dataset[i].measSeq); } FREE (dataset); @@ -90,14 +94,75 @@ } +static int NsynthNames = 6; +static float synthMaxMags[] = {13.5, 13.5, 13.5, 13.0, 12.0, 13.5}; +static char *synthNames[] = { "g", "r", "i", "z", "y", "w"}; + +void SetMrelInfoResetObject (SetMrelInfo *results) { + + int i; + for (i = 0; i < results->Nsecfilt; i++) { + results->havePS1[i] = FALSE; + results->haveSYN[i] = FALSE; + results->needSYN[i] = FALSE; + results->measSYN[i] = -1; + results->minSYN[i] = NAN; + + results->psfQfMax[i] = NAN; + results->psfQfPerfMax[i] = NAN; + + results->Next[i] = 0; + results->NexpPS1[i] = 0; + results->havePS1[i] = 0; + results->haveUbercal[i] = 0; + results->minUbercalDist[i] = 1000; + + results->psfData[i].Nlist = 0; + results->aperData[i].Nlist = 0; + results->kronData[i].Nlist = 0; + } + results->NextPS1 = 0; + results->NpsfPS1 = 0; +} + void SetMrelInfoInit (SetMrelInfo *results, int allocLists) { + + int i; SetMrelInfoReset (results); if (allocLists) { results->Nsecfilt = GetPhotcodeNsecfilt (); + results->psfData = StatDataSetAlloc (results->Nsecfilt, Nmax); results->aperData = StatDataSetAlloc (results->Nsecfilt, Nmax); results->kronData = StatDataSetAlloc (results->Nsecfilt, Nmax); + ALLOCATE (results->psfqf_list, double, Nmax); ALLOCATE (results->psfqfperf_list, double, Nmax); ALLOCATE (results->stargal_list, double, Nmax); + + ALLOCATE (results->havePS1, int, results->Nsecfilt); + ALLOCATE (results->haveSYN, int, results->Nsecfilt); + ALLOCATE (results->needSYN, int, results->Nsecfilt); + ALLOCATE (results->measSYN, int, results->Nsecfilt); + ALLOCATE (results->minSYN, float, results->Nsecfilt); + + ALLOCATE (results->psfQfMax, float, results->Nsecfilt); + ALLOCATE (results->psfQfPerfMax, float, results->Nsecfilt); + + ALLOCATE (results->Next, int, results->Nsecfilt); + ALLOCATE (results->NexpPS1, int, results->Nsecfilt); + ALLOCATE (results->havePS1, int, results->Nsecfilt); + ALLOCATE (results->haveUbercal, int, results->Nsecfilt); + + ALLOCATE (results->minUbercalDist, float, results->Nsecfilt); + + SetMrelInfoResetObject (results, results->Nsecfilt); + + for (i = 0; i < NsynthNames; i++) { + int photcode = GetPhotcodeCodebyName (synthNames[i]); + if (!photcode) continue; + int Nsec = GetPhotcodeNsec (photcode); + if (Nsec < 0) continue; + results->minSYN[Nsec] = synthMaxMags[i]; + } } else { results->Nsecfilt = 0; @@ -105,17 +170,32 @@ results->aperData = NULL; results->kronData = NULL; + results->psfqf_list = NULL; results->psfqfperf_list = NULL; results->stargal_list = NULL; + + results->havePS1 = NULL; + results->haveSYN = NULL; + results->needSYN = NULL; + results->measSYN = NULL; + results->minSYN = NULL; + + results->psfQfMax = NULL; + results->psfQfPerfMax = NULL; + results->Next = NULL; + results->NexpPS1 = NULL; + results->havePS1 = NULL; + results->haveUbercal = NULL; + results->minUbercalDist = NULL; } } void SetMrelInfoFree (SetMrelInfo *results) { - StatDataSetFree (results->psfdata, results->Nsecfilt); - StatDataSetFree (results->aperdata, results->Nsecfilt); - StatDataSetFree (results->krondata, results->Nsecfilt); - FREE (results.psfqf_list); - FREE (results.psfqfperf_list); - FREE (results.stargal_list); + StatDataSetFree (results->psfData, results->Nsecfilt); + StatDataSetFree (results->aperData, results->Nsecfilt); + StatDataSetFree (results->kronData, results->Nsecfilt); + FREE (results->psfqf_list); + FREE (results->psfqfperf_list); + FREE (results->stargal_list); } Index: /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/extra.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/extra.c (revision 38292) +++ /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/extra.c (revision 38293) @@ -70,2 +70,9 @@ return FALSE; } + +int is2MASS (int photcode) { + + if ((photcode >= 2011) && (photcode <= 2013)) return TRUE; + return FALSE; +} + Index: /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/liststats.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/liststats.c (revision 38292) +++ /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/liststats.c (revision 38293) @@ -14,4 +14,19 @@ fprintf (stderr, "ERROR: invalid stats mode: %s\n", strmode); exit (1); +} + +int liststats_init (StatType *stats) { + stats->median = NAN; + stats->mean = NAN; + stats->sigma = NAN; + stats->error = NAN; + stats->chisq = NAN; + stats->min = NAN; + stats->max = NAN; + stats->Upper80 = NAN; + stats->Lower20 = NAN; + stats->total = NAN; + stats->Nmeas = 0; + return TRUE; } Index: /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelCatalog.c (revision 38292) +++ /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelCatalog.c (revision 38293) @@ -2,10 +2,10 @@ # define UBERCAL_WEIGHT 100.0 -# define SKIP_THIS_MEAS(REASON) { \ - measureT[k].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ +# define SKIP_THIS_MEAS(REASON) { \ + measureT[k].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ if (measure) { \ - measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ - } \ - results->REASON ++; \ + measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ + } \ + results->REASON ++; \ continue; } @@ -13,4 +13,6 @@ measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ continue; } + +# define CHECK_VALID_MAG(MAG) (!isnan(MAG) && (MAG > 0.0) && (MAG < 30.0)) static float MagToFlux (float Mag) { @@ -54,5 +56,5 @@ Measure *measure = catalog[Nc].measure ? &catalog[Nc].measure[m] : NULL; SecFilt *secfilt = catalog[Nc].secfilt ? &catalog[Nc].secfilt[j*Nsecfilt] : NULL; - setMrelAverageExposure (m, Nc, pass, flatcorr, results, average, &catalog[Nc].averageT[j], secfilt, measure, &catalog[Nc].measureT[m], &catalog[Nc].found_t[Nsecfilt*j]); + setMrelAverageExposure (m, Nc, isSetMrelFinal, flatcorr, results, average, &catalog[Nc].averageT[j], secfilt, measure, &catalog[Nc].measureT[m], &catalog[Nc].found_t[Nsecfilt*j]); // only apply Stack operation on setMrelFinal in first pass @@ -68,10 +70,314 @@ } } - // if (primaryCell) free (primaryCell); + return (TRUE); } +// NOTE: +// Msys is measure[i].M + zp corrections +// Mcal is image[j].Mcal +// Mmos and Mgrid are offsets for mosaic and grid + +// 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 +// XXX probably need to use the photocde table to assign reference mag weights. + +// 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 + // set mean of chip measurements (selected by photcode range for now): -int setMrelAverageExposure (off_t measureOffset, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found) { +int setMrelAverageExposure (off_t measureOffset, int cat, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT) { + + off_t k; + + // we are measuring means for 3 types of mags: psf, ap, kron. + + // in the final assignment, set the mean mags even if only 1 measurements exists + int Nminmeas = isSetMrelFinal ? 1 : STAR_TOOFEW + 1; + + // isSetMrelFinal : in the final pass, we set psf, kron, ap mag values and extra stats; other passes only do psf mags + + // option for a test print + if (FALSE && (averageT[0].objID == 0x7146) && (averageT[0].catID == 0x49d8)) { + fprintf (stderr, "test obj\n"); + print_measure_set_alt (average, secfilt, measure); + } + + SetMrelInfoResetObject (results); // reset the per-object arrays + + // assign measurements to the photcode lists + for (k = 0; k < averageT[0].Nmeasure; k++) { + + // skip measurements that do not have a valid photcode (raise exception?) + PhotCode *code = GetPhotcodebyCode (measureT[k].photcode); + myAssert (code, "invalid photcode??"); + + // SKIP gpc1 stack data (hard-wired photcodes) + if (isGPC1stack(measureT[k].photcode)) continue; + + // SKIP gpc1 forced-warp data (hard-wired photcodes) + if (isGPC1warp(measureT[k].photcode)) continue; + + if (isTycho(measureT[k].photcode)) { haveTYCHO = TRUE; } + + if (is2MASS(measureT[k].photcode)) { + if (measureT[k].photFlags & 0x00c00000) { + Galaxy2MASS = TRUE; // per object value + } + if (measureT[k].photFlags & 0x00000007) { + Good2MASS = TRUE; // per object value + } + } + + // assign the Nsec value so we can assign to the right lists + int Nsec = GetPhotcodeNsec (code->equiv); + if (Nsec < 0) continue; // skip measurements which do not have an equiv average photcode + + // force the use of SYN even if we have PS1 mags? + if (isGPC1synth(measureT[k].photcode)) { + results->haveSYN[Nsec] = TRUE; + results->measSYN[Nsec] = k; + if (measureT[k].M < results->minSYN[Nsec]) results->needSYN[Nsec] = TRUE; + continue; + } + + // assign the max psfQF values for each Nsec value + if (measureT[k].psfQF > results->psfQfMax[Nsec]) results->psfQfMax[Nsec] = measureT[k].psfQF; + if (measureT[k].psfQFperf > results->psfQfPerfMax[Nsec]) results->psfQfPerfMax[Nsec] = measureT[k].psfQFperf; + + // are we a PS1 exposure photcode? (hard-wired photcodes) + if (isGPC1chip(measureT[k].photcode)) { results->NexpPS1[Nsec] ++; results->havePS1[Nsec] = TRUE; } + + // ** Choose the calibration (depends on the mode : do I have an image reference or not?) + float Mcal = 0, Mmos = 0, Mgrid = 0; + if (getImageEntry (measureOffset + k, cat) < 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 = measureT[k].Mcal; // check that this is zero for loaded REF value + } else { + Mcal = getMcal_alt (meas, cat, flatcorr, measureT[k].Xccd, measureT[k].Yccd); + if (isnan(Mcal)) SKIP_THIS_MEAS(Ncal); + Mmos = getMmos (meas, cat); + if (isnan(Mmos)) SKIP_THIS_MEAS(Nmos); + Mgrid = getMgrid (meas, cat); + if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid); + } + + // XXX ubercal dist should be by Nsec + int myUbercalDist = getUbercalDist(meas, cat); + result->minUbercalDist[Nsec] = MIN(result->minUbercalDist[Nsec], myUbercalDist); + + int isUbercal = (measureT[k].dbFlags & ID_MEAS_PHOTOM_UBERCAL); + + if (isUbercal) results->haveUbercal[Nsec] = TRUE; // haveUbercal is set per secfilt, isUbercal is per measure XXX define this array + + int useUbercalWeight = isUbercal || (refPhotcode && (code->code == refPhotcode->code)); + + float Map, dMap; + if (isSetMrelFinal) { + Map = PhotCat (&measure[k], MAG_CLASS_APER); + dMap = MAX (hypot(measureT[k].dMap, code->photomErrSys), MIN_ERROR); + if (CHECK_VALID_MAG(Map)) { + // XXX need to carry a counter for each list set + int Nap = results->apData[Nsec].Nlist; + results->apData[Nsec].flxlist[Nap] = Map - Mcal - Mmos - Mgrid; + results->apData[Nsec].errlist[Nap] = dMpsf; // no other option + results->apData[Nsec].wgtlist[Nap] = useUbercalWeight ? UBERCAL_WEIGHT : 1.0; + results->apData[Nsec].ranking[Nap] = measureRank[k]; + results->apData[Nsec].measSeq[Nap] = k; + results->apData[Nsec].Nlist ++; + } + } + + float Mkron, dMkron; + if (isSetMrelFinal) { + Mkron = PhotCat (&measure[k], MAG_CLASS_KRON); + dMkron = MAX (hypot(measureT[k].dMkron, code->photomErrSys), MIN_ERROR); + if (CHECK_VALID_MAG(Mkron)) { + int Nkron = results->kronData[Nsec].Nlist; + results->kronData[Nsec].flxlist[Nkron] = Mkron - Mcal - Mmos - Mgrid; + results->kronData[Nsec].errlist[Nkron] = measure[k].dMkron; + results->kronData[Nsec].wgtlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0; + results->kronData[Nsec].ranking[Nkron] = measureRank[k]; + results->kronData[Nsec].measSeq[Nkron] = k; + results->kronData[Nsec].Nlist ++; + } + } + + float Mpsf = PhotSysTiny (&measureT[k], &averageT[0], &secfilt[0], MAG_CLASS_PSF); + float dMpsf = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR); + if (CHECK_VALID_MAG(Mpsf)) { + int Npsf = results->psfData[Nsec].Nlist; + results->psfData[Nsec].flxlist[Npsf] = Mpsf - Mcal - Mmos - Mgrid; + results->psfData[Nsec].errlist[Npsf] = dMpsf; + results->psfData[Nsec].wgtlist[Npsf] = isUbercal ? UBERCAL_WEIGHT : 1.0; + results->psfData[Nsec].ranking[Npsf] = measureRank[k]; + results->psfData[Nsec].measSeq[Npsf] = k; + results->psfData[Nsec].Nlist ++; + } else { + measureT[k].dbFlags |= ID_MEAS_SKIP_PHOTOM; // XXX is this used outside of relphot? + if (measure) measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM; + results->Nsys ++; + } + + // count the extended detections (all PS1 bands) + if (isGPC1chip(measureT[k].photcode) && !isnan(Map) && !isnan(Mpsf)) { + float dMagAp = Mpsf - Map; + float SigmaAp = hypot(0.1, 2.5*dMpsf); + // 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) { + results->Next[Nsec] ++; + results->NextPS1 ++; + } else { + results->NpsfPS1 ++; + } + } + } + + // now calculate the mean stats for the Nsec bands. + for (Nsec = 0; Nsec < Nsecfilt; Nsec++) { + dvo_init_secfilt (&secfilt[Nsec]); + + // XXX hardwired + if ((thisCode < 6) || (thisCode == 9)) { + secfilt[Nsec].Ncode = results->NexpPS1[Nsec]; + } else { + secfilt[Nsec].Ncode = results->Ncode[Nsec]; // 2MASS data if it exists + } + + // XXX need to also supply ZP from map, apply tycho logic + if (results->needSYN[Nsec]) { + if (!isSetMrelFinal) continue; + + // use the single SYNTH value instead of the other mags here + float Mpsf = measure[results->measSYN[Nsec]].M; + float ZP = 0.0; + + if (!haveTycho) { + // need to look up the (X,Y) coords from this (R,D) location + double X, Y; + double R = ohana_normalize_angle_to_midpoint (measure->R, 0.0); // XXX 0.0 or 180.0 for center? + RD_to_XY (&X, &Y, R, measure->D, &zpts->coords); + if (X < 0) return FALSE; + if (Y < 0) return FALSE; + if (X >= zpts->Nx) return FALSE; + if (Y >= zpts->Ny) return FALSE; + + int Xpix = X; + int Ypix = Y; + int Npix = Xpix + Ypix*zpts->Nx; + + float *value = (float *) zpts->matrix[Nsec].buffer; + ZP = !isnan(value[Npix]) ? value[Npix] : 0.0; + secfilt[Nsec].flags |= ID_SECF_HAVE_TYCHO; + } + secfilt[Nsec].M = Mpsf + ZP; + secfilt[Nsec].dM = 0.6; + secfilt[Nsec].chipsq = 0.0; + secfilt[Nsec].flags |= ID_SECF_USE_SYNTH; + continue; + } + + // if too few valid measurements meet the minimum criteria, go to the next entry + int Nranking = magStatsByRanking (&results->psfData[Nsec], &results->psfstats); + if (Nranking < Nminmeas) { + secfilt[Nsec].flags |= ID_STAR_FEW; + } else { + secfilt[Nsec].M = psfstats->mean; + secfilt[Nsec].dM = psfstats->error; + secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN; + } + int minRank = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10; + + // when running -averages, we have no information about the images, so we cannot set this + if (minUbercalDist > -1) { + secfilt[Nsec].ubercalDist = minUbercalDist; + } + + if (isSetMrelFinal) { + secfilt[Nsec].Mstdev = psfstats->sigma; + secfilt[Nsec].Nused = psfstats->Nmeas; + secfilt[Nsec].Mmax = psfstats->max; + secfilt[Nsec].Mmin = psfstats->min; + + secfilt[Nsec].psfQfMax = psfQfMax; + secfilt[Nsec].psfQfPerfMax = psfQfPerfMax; + + Nranking = magStatsByRanking (&results->aperData[Nsec], results->apstats); + secfilt[Nsec].Map = apstats->mean; + secfilt[Nsec].dMap = apstats->error; + secfilt[Nsec].sMap = apstats->sigma; + secfilt[Nsec].NusedAp = Nranking; + + Nranking = magStatsByRanking (&results->kronData[Nsec], kronstats); + secfilt[Nsec].Mkron = kronstats->mean; + secfilt[Nsec].dMkron = kronstats->error; + secfilt[Nsec].sMkron = kronstats->sigma; + secfilt[Nsec].NusedKron = Nranking; + + // does this object appear extended in > 50% of measurements? + if ((results->Next[Nsec] > 0) && (results->Next[Nsec] > 0.5*results->Npsf[Nsec])) { + secfilt[Nsec].flags |= ID_SECF_OBJ_EXT; + } + + switch (minRank) { + case 0: + secfilt[Nsec].flags |= ID_PHOTOM_PASS_0; + if (results->havePS1[Nsec]) GoodPS1 = TRUE; + break; + case 1: + secfilt[Nsec].flags |= ID_PHOTOM_PASS_1; + if (results->havePS1[Nsec]) GoodPS1 = TRUE; + break; + case 2: + secfilt[Nsec].flags |= ID_PHOTOM_PASS_2; + if (results->havePS1[Nsec]) GoodPS1 = TRUE; + break; + case 3: + secfilt[Nsec].flags |= ID_PHOTOM_PASS_3; + break; + case 4: + secfilt[Nsec].flags |= ID_PHOTOM_PASS_4; + break; + default: + break; + } + if (results->havePS1[Nsec]) { + secfilt[Nsec].flags |= ID_SECF_HAS_PS1; + } + if (results->haveUbercal[Nsec]) { + secfilt[Nsec].flags |= ID_SECF_USE_UBERCAL; + } + } + } + + if (isSetMrelFinal) { + if (results->NextPS1 && (results->NextPS1 > results->NpsfPS1)) { + average[0].flags |= ID_OBJ_EXT; + } + if (GoodPS1) { + average[0].flags |= ID_OBJ_GOOD; + } + if (Galaxy2MASS) { + average[0].flags |= ID_OBJ_EXT_ALT; + } + if (Good2MASS) { + average[0].flags |= ID_OBJ_GOOD_ALT; + } + } + return (TRUE); +} + +// set mean of chip measurements (selected by photcode range for now): +int setMrelAverageExposure_old (off_t measureOffset, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found) { off_t k; @@ -81,17 +387,15 @@ // error for the ap mags, but krons have their own errors. it is an open question if I // should be doing weighted or unweighted fits (this is a user option) - - // these lists are pre-allocated to the max length - double **Mpsflist = results->Mpsflist; - double **dpsflist = results->dpsflist; - double **wpsflist = results->wpsflist; - - double **Maplist = results->Maplist; - double **daplist = results->daplist; - double **waplist = results->waplist; - - double **Mkronlist = results->Mkronlist; - double **dkronlist = results->dkronlist; - double **wkronlist = results->wkronlist; + double *Mpsflist = results->Mpsflist; + double *dpsflist = results->dpsflist; + double *wpsflist = results->wpsflist; + + double *Maplist = results->Maplist; + double *daplist = results->daplist; + double *waplist = results->waplist; + + double *Mkronlist = results->Mkronlist; + double *dkronlist = results->dkronlist; + double *wkronlist = results->wkronlist; StatType *psfstats = &results->psfstats; @@ -105,62 +409,4 @@ fprintf (stderr, "test obj\n"); print_measure_set_alt (average, secfilt, measure); - } - - // first pass: check for PS1[n], Synth, Tycho, 2MASS - - // first, determine the state of things for each of the main filters: - // in this bit, we are using hardwired photcode values: - // PS1_g, PS1_r, PS1_i, PS1_z, PS1_y, PS1_w == (1,2,3,4,5,6); 0 == none - // note that Nsec = PS_N - 1 - - int havePS1[7]; - int haveSYN[7]; - int needSYN[7]; - int measSYN[7]; - - int haveTYCHO = FALSE; - int have2MASS = FALSE; - - // was: 13.64, 13.76, 13.74, 12.94, 12.01 - float MaxMagForceSynth[] = {0.0, 13.5, 13.5, 13.5, 13.0, 12.0, 13.5}; - // NONE, g, r, i, z, y, w - - for (i = 0; i < 7; i++) { - havePS1[i] = FALSE; // do we have any PS1 measurement for this filter? - haveSYN[i] = FALSE; // do we have any synthetic measurement for this filter? - needSYN[i] = FALSE; // force the use of the synthetic measurement for this filter? - measSYN[i] = -1; // which entries carries the synthetic measurement for this filter? - } - - for (k = 0; k < average[0].Nmeasure; k++) { - - // skip measurements that do not have a valid photcode (which ones lack photcode?) - PhotCode *code = GetPhotcodebyCode (measure[k].photcode); - if (!code) continue; - - // is this synth_grizy? - // is this PS1? - - // are we a PS1 exposure photcode? - if (isGPC1chip(measure[k].photcode)) { - int Nfilter = whichGPC1filter (measure[k].photcode); - havePS1[Nfilter] = TRUE; - } - if (isGPC1synth(measure[k].photcode)) { - int Nfilter = measure[k].photcode - 3000; - haveSYN[Nfilter] = TRUE; - measSYN[Nfilter] = k; - - // force the use of SYN even if we have PS1 mags? - if (measure[k].M < MaxMagForceSynth[Nfilter]) { - needSYN[Nfilter] = TRUE; - } - } - if (isTycho(measure[k].photcode)) { - haveTYCHO = TRUE; - } - if (is2MASS(measure[k].photcode)) { - have2MASS = TRUE; - } } @@ -505,377 +751,4 @@ } -// set mean of chip measurements (selected by photcode range for now): -int setMrelAverageExposure_old (off_t measureOffset, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found) { - - off_t k; - float Msys = 0, Mcal= 0, Mmos = 0, Mgrid = 0; - - // we are measuring means for 3 types of mags: psf, ap, kron. I am using the psf mag - // error for the ap mags, but krons have their own errors. it is an open question if I - // should be doing weighted or unweighted fits (this is a user option) - double *Mpsflist = results->Mpsflist; - double *dpsflist = results->dpsflist; - double *wpsflist = results->wpsflist; - - double *Maplist = results->Maplist; - double *daplist = results->daplist; - double *waplist = results->waplist; - - double *Mkronlist = results->Mkronlist; - double *dkronlist = results->dkronlist; - double *wkronlist = results->wkronlist; - - StatType *psfstats = &results->psfstats; - StatType *apstats = &results->apstats; - StatType *kronstats = &results->kronstats; - - int isSetMrelFinal = (pass >= 0); - - // option for a test print - if (FALSE && (average[0].objID == 0x7146) && (average[0].catID == 0x49d8)) { - fprintf (stderr, "test obj\n"); - print_measure_set_alt (average, secfilt, measure); - } - - int GoodPS1 = FALSE; - int Good2MASS = FALSE; - int Galaxy2MASS = FALSE; - - int NextPS1 = 0; - int NpsfPS1 = 0; - - 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 && found[Nsec]) continue; - - // skip bad stars - if (!isSetMrelFinal && (secfilt[Nsec].flags & STAR_BAD)) continue; - - int NexpPS1 = 0; - int Ncode = 0; - int Next = 0; - int haveSynth = FALSE; - int havePS1 = FALSE; - - int forceSynth = FALSE; - int forceSynthEntry = -1; - - int haveUbercal = FALSE; - - int minUbercalDist = 1000; - - off_t meas = measureOffset; - - float psfQfMax = 0.0; - float psfQfPerfMax = 0.0; - - int Nap = 0; - int Npsf = 0; - int Nkron = 0; - for (k = 0; k < averageT[0].Nmeasure; k++, meas++) { - - // skip measurements that do not match the current photcode - PhotCode *code = GetPhotcodebyCode (measureT[k].photcode); - if (!code) continue; - if (code->equiv != thisCode) { continue; } - Ncode ++; - - // are we a PS1 exposure photcode? - if (isGPC1chip(measureT[k].photcode)) NexpPS1 ++; - - // SKIP gpc1 stack data - if (isGPC1stack(measureT[k].photcode)) continue; - - // SKIP gpc1 forced-warp data - if (isGPC1warp(measureT[k].photcode)) continue; - - if (isSetMrelFinal) { - if (measure[k].psfQF > psfQfMax) psfQfMax = measure[k].psfQF; - if (measure[k].psfQFperf > psfQfPerfMax) psfQfPerfMax = measure[k].psfQFperf; - } - - if (measureT[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); - - if (getImageEntry (meas, cat) < 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 = measureT[k].Mcal; // check that this is zero for loaded REF value - } else { - Mcal = getMcal_alt (meas, cat, flatcorr, measureT[k].Xccd, measureT[k].Yccd); - if (isnan(Mcal)) SKIP_THIS_MEAS(Ncal); - Mmos = getMmos (meas, cat); - if (isnan(Mmos)) SKIP_THIS_MEAS(Nmos); - Mgrid = getMgrid (meas, cat); - if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid); - } - - // skip some absurd values NAN, < 0.0, > 30.0 - Msys = PhotSysTiny (&measureT[k], &averageT[0], &secfilt[0], 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(meas, cat); - minUbercalDist = MIN(minUbercalDist, myUbercalDist); - - int isUbercal = (measureT[k].dbFlags & ID_MEAS_PHOTOM_UBERCAL); - float dMpsf = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR); - - if (isUbercal) haveUbercal = TRUE; // haveUbercal is set per secfilt, isUbercal is per measure - - if (isSetMrelFinal) { - // special options for PS1 data - if ((measure[k].photcode >= 10000) && (measure[k].photcode <= 10500)) { - // count the extended detections - if (!isnan(measure[k].Map)) { - float dMagAp = measure[k].M - measure[k].Map; - float SigmaAp = hypot(0.1, 2.5*measure[k].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; - } - - // count extended detections for 2MASS (XXX NOTE hardwired photcodes 2011, 2012, 2013) - if ((measure[k].photcode >= 2011) && (measure[k].photcode <= 2013)) { - if (measure[k].photFlags & 0x00c00000) { - Next ++; - Galaxy2MASS = TRUE; - } - if (pass == 0) { - if (measure[k].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 ((measure[k].photcode >= 3001) && (measure[k].photcode <= 3005)) { - // something of a hack: force object to use synth values if synth mags >> - // saturation (3pi instrumental mags < -15) - float MaxMagForceSynth = NAN; - switch (measure[k].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 (measureT[k].M < MaxMagForceSynth) { - forceSynth = TRUE; - forceSynthEntry = Npsf; - } else { - if (pass < 4) { - SKIP_THIS_MEAS(Nbad); - } - haveSynth = TRUE; - } - } - - // Map (Maplist) and Mkron (Mkronlist,dkronlist) are used to calculate mean mags per filter - float Map = PhotCat (&measure[k], MAG_CLASS_APER); - if (!isnan(Map)) { - Maplist[Nap] = Map - Mcal - Mmos - Mgrid; - daplist[Nap] = dMpsf;// XXX check on this... - waplist[Nap] = isUbercal ? UBERCAL_WEIGHT : 1.0; - Nap ++; - } - - float Mkron = PhotCat (&measure[k], MAG_CLASS_KRON); - if (!isnan(Mkron)) { - Mkronlist[Nkron] = Mkron - Mcal - Mmos - Mgrid; - dkronlist[Nkron] = measure[k].dMkron; - wkronlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0; - Nkron ++; - } - } // if (isSetMrelFinal) - - // 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 - Mpsflist[Npsf] = Msys - Mcal - Mmos - Mgrid; - dpsflist[Npsf] = dMpsf; - wpsflist[Npsf] = isUbercal ? UBERCAL_WEIGHT : 1.0; - - // NOTE: - // Msys is measure[i].M + zp corrections - // Mcal is image[j].Mcal - // Mmos and Mgrid are offsets for mosaic and grid - - // 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 - // XXX probably need to use the photocde table to assign reference mag weights. - if (refPhotcode) { - if (code->code == refPhotcode->code) { - wpsflist[Npsf] = UBERCAL_WEIGHT; - } - } - Npsf ++; - } - - int Nminmeas = isSetMrelFinal ? 1 : STAR_TOOFEW + 1; - - // XXX : ugh : another hard-wired photcode entry... - if (isSetMrelFinal && (pass == 0)) { - if ((thisCode < 6) || (thisCode == 9)) { - secfilt[Nsec].Ncode = NexpPS1; - } else { - secfilt[Nsec].Ncode = Ncode; // 2MASS data if it exists - } - } - - // when performing the grid analysis, STAR_TOOFEW should be set to 1; - if (Npsf < Nminmeas) { /* too few measurements */ - // fprintf (f, "%10.6f %10.6f %d %d %d\n", averageT[0].R, averageT[0].D, measureT[0].imageID, Npsf, STAR_TOOFEW); - secfilt[Nsec].flags |= ID_STAR_FEW; - if (Ncode == 0) { - results->Ncode ++; - } else { - results->Nfew ++; - } - continue; - } else { - secfilt[Nsec].flags &= ~ID_STAR_FEW; - } - - if (forceSynth) { - // use the single SYNTH value instead of the other mags here - myAssert ((forceSynthEntry < Npsf) && (forceSynthEntry >= 0), "programming error"); - Mpsflist[0] = Mpsflist[forceSynthEntry]; - dpsflist[0] = dpsflist[forceSynthEntry]; - wpsflist[0] = wpsflist[forceSynthEntry]; - Npsf = 1; - } - liststats (Mpsflist, dpsflist, wpsflist, Npsf, psfstats); - - secfilt[Nsec].M = psfstats->mean; - secfilt[Nsec].dM = psfstats->error; - secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN; - - // when running -averages, we have no information about the images, so we cannot set this - if (minUbercalDist > -1) { - secfilt[Nsec].ubercalDist = minUbercalDist; - } - - if (isSetMrelFinal) { - found[Nsec] = TRUE; - - secfilt[Nsec].Mstdev = psfstats->sigma; // Mstdev is in millimags (not enough space for more precision) - // secfilt[Nsec].Ncode = Ncode; - secfilt[Nsec].Nused = psfstats->Nmeas; - - secfilt[Nsec].Mmax = psfstats->max; - secfilt[Nsec].Mmin = psfstats->min; - - secfilt[Nsec].psfQfMax = psfQfMax; - secfilt[Nsec].psfQfPerfMax = psfQfPerfMax; - - // NOTE : use the modified weight for apmags as well as psf mags - liststats (Maplist, daplist, waplist, Nap, apstats); - secfilt[Nsec].Map = Nap > 0 ? apstats->mean : NAN; - secfilt[Nsec].dMap = Nap > 0 ? apstats->error : NAN; - secfilt[Nsec].sMap = Nap > 0 ? apstats->sigma : NAN; - secfilt[Nsec].NusedAp = Nap; - - liststats (Mkronlist, dkronlist, wkronlist, Nkron, kronstats); - secfilt[Nsec].Mkron = Nkron > 0 ? kronstats->mean : NAN; - secfilt[Nsec].dMkron = Nkron > 0 ? kronstats->error : NAN; - secfilt[Nsec].sMkron = Nkron > 0 ? kronstats->sigma : NAN; - secfilt[Nsec].NusedKron = Nkron; - - // NOTE: for 2MASS measurements, Next should be 1, as should N - if ((Next > 0) && (Next > 0.5*Npsf)) { - secfilt[Nsec].flags |= ID_SECF_OBJ_EXT; - } - - switch (pass) { - case 0: - secfilt[Nsec].flags |= ID_PHOTOM_PASS_0; - if (havePS1) GoodPS1 = TRUE; - break; - case 1: - secfilt[Nsec].flags |= ID_PHOTOM_PASS_1; - if (havePS1) GoodPS1 = TRUE; - break; - case 2: - secfilt[Nsec].flags |= ID_PHOTOM_PASS_2; - if (havePS1) GoodPS1 = TRUE; - break; - case 3: - secfilt[Nsec].flags |= ID_PHOTOM_PASS_3; - break; - case 4: - secfilt[Nsec].flags |= ID_PHOTOM_PASS_4; - break; - } - if (haveSynth) { - secfilt[Nsec].flags |= ID_SECF_USE_SYNTH; - } - if (havePS1) { - secfilt[Nsec].flags |= ID_SECF_HAS_PS1; - } - if (haveUbercal) { - secfilt[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: - average[0].flags &= ~flagBits; - } - - if (NextPS1 && (NextPS1 > NpsfPS1)) { - average[0].flags |= ID_OBJ_EXT; - } - if (GoodPS1) { - average[0].flags |= ID_OBJ_GOOD; - } - if (Galaxy2MASS) { - average[0].flags |= ID_OBJ_EXT_ALT; - } - if (Good2MASS) { - average[0].flags |= ID_OBJ_GOOD_ALT; - } - } - return (TRUE); -} - // only apply Stack operation on setMrelFinal in first pass // this function has 3 goals, not to be confused: @@ -1338,2 +1211,28 @@ return (TRUE); } + +int magStatsByRanking (StatDataSet *dataset, StatType *stats) { + + liststats_init (stats); + + if (dataset->Nlist == 0) return 0; + + // we have a list of measurements for this Nsec value, along with errors, weights, and their rank + // first, sort by the rank (increasing) + sort_measures_by_ranking (dataset->ranking, dataset->measSeq, dataset->flxlist, dataset->errlist, dataset->wgtlist, dataset->Nlist); + + // find the values with the same minimum rank: + int Nranking = 0; + int minRank = dataset->ranking[0]; // MIN (5, result->psfData[Nsec].ranking[0]); -- only allow rank 5 or < + + int i; + for (i = 0; (i < dataset->Nlist) && (dataset->ranking[i] == minRank); i++, Nranking++); + + // if no valid measurements meet the minimum criteria, go to the next entry + if (Nranking < Nminmeas) { + return 0; + } + + liststats (dataset->flxlist, dataset->errlist, dataset->wgtlist, Nranking, stats); + return (Nranking); +} Index: /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelFinal.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelFinal.c (revision 38292) +++ /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelFinal.c (revision 38293) @@ -12,5 +12,5 @@ int Nsecfilt = GetPhotcodeNsecfilt (); - /* if we reset the catalog, reset all the current measurements */ + /*** RESET photometry flags and Mcal values ***/ if (RESET) { // flags used by the photometry analysis (excluding UBERCAL) @@ -100,8 +100,8 @@ // XXX make this optional? (do not clean for -averages?) if (!simpleAverage) clean_measures (catalog, 1, TRUE, flatcorr); /* mark outliers ID_MEAS_POOR_PHOTOM */ - for (i = 0; i < 5; i++) { - skip_measurements (catalog, i, flatcorr); /* set ID_MEAS_SKIP for measures to be skipped */ - setMrelOutput (catalog, 1, i, flatcorr); /* set Mrel using remaining measures */ - } + + ALLOCATE (catalog[0].measureRank, char, catalog[0].Nmeasure); + setMeasurementRank (catalog); + setMrelOutput (catalog); if (!simpleAverage) setMcalOutput (catalog, 1, flatcorr); @@ -267,14 +267,79 @@ } -/* dophot types: - -1 - star -2 - galaxy -3 - star -4 - weak -5 - edge -7 - weak -9 - weak -10 - sat -*/ - +void setMeasureRank (Catalog *catalog) { + + int i; + + // make these global + int IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL; + int IMAGE_POOR = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP; + int MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM; + int MEAS_POOR = ID_MEAS_POOR_PHOTOM | ID_MEAS_AREA; + + Measure *measure = catalog[0].measure; + MeasureTiny *measureT = catalog[0].measureT; + char *measureRank = catalog[0].measureRank; + + /* set measureRank[] based on various quality measurements */ + for (i = 0; i < catalog[0].Nmeasure; i++) { + measureRank[i] = 11; // start at a low rank + + // clear SKIP for all measures at first + measureT[i].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; + measure [i].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; + + // skip measurements without a valid photcode + PhotCode *code = GetPhotcodebyCode (measure[i].photcode); + if (!code) continue; + + // measurements outside time range have poor rank + if (TimeSelect) { + if (measure[i].t < TSTART) { measureRank[i] = 10; continue; } + if (measure[i].t > TSTOP) { measureRank[i] = 10; continue; } + } + + // measurements from BAD images and mosaics (not REF mags) + off_t Nim = getImageEntry (i, 0); + int imageFlags = getImageFlags (i, 0); + int mosaicFlags = getMosaicFlags (i, 0); + + if (Nim > -1) { + if (imageFlags & IMAGE_BAD) { measureRank[i] = 9; continue; } + + if (mosaicFlags & IMAGE_BAD) { measureRank[i] = 9; continue; } + + // measurements ranked by inst mag limit (not REF mags) + if (ImagSelect) { + mag = PhotInst (&measure[i], MAG_CLASS_PSF); + if (mag < ImagMin) { measureRank[i] = 8; continue; } + if (mag > ImagMax) { measureRank[i] = 8; continue; } + } + } + + // RANK 1 : BAD photFlags + if (measure[i].dbFlags & MEAS_BAD) { measureRank[i] = 7; continue; } + + // RANK 3 : BAD photFlags (eg, SAT, CR), internal outliers + if (measure[i].photFlags & code->photomBadMask) { measureRank[i] = 6; continue; } + + if (Nim > -1) { + if (imageFlags & IMAGE_POOR) { measureRank[i] = 5; continue; } + if (mosaicFlags & IMAGE_POOR) { measureRank[i] = 5; continue; } + } + + // RANK 2 : psfQF value + if (measure[i].psfQF < 0.85) { measureRank[i] = 4; continue; } + + // RANK 1 : POOR photFlags + if (measure[i].dbFlags & MEAS_POOR) { measureRank[i] = 3; continue; } + + // RANK 1 : POOR photFlags + if (measure[i].photFlags & code->photomPoorMask) { measureRank[i] = 2; continue; } + + // RANK 1 : psfQFperfect value + if (measure[i].psfQFperf < 0.85) { measureRank[i] = 1; continue; } + + // RANK 0 : perfect measurement: + measureRank[i] = 0; + } +}