Index: /branches/eam_branches/ipp-20121130/Ohana/src/relphot/Makefile =================================================================== --- /branches/eam_branches/ipp-20121130/Ohana/src/relphot/Makefile (revision 34836) +++ /branches/eam_branches/ipp-20121130/Ohana/src/relphot/Makefile (revision 34837) @@ -49,4 +49,5 @@ $(SRC)/setExclusions.$(ARCH).o \ $(SRC)/setMrelFinal.$(ARCH).o \ +$(SRC)/setMrelCatalog.$(ARCH).o \ $(SRC)/BoundaryTreeOps.$(ARCH).o \ $(SRC)/write_coords.$(ARCH).o @@ -78,4 +79,5 @@ $(SRC)/setExclusions.$(ARCH).o \ $(SRC)/setMrelFinal.$(ARCH).o \ +$(SRC)/setMrelCatalog.$(ARCH).o \ $(SRC)/BoundaryTreeOps.$(ARCH).o \ $(SRC)/write_coords.$(ARCH).o Index: /branches/eam_branches/ipp-20121130/Ohana/src/relphot/include/relphot.h =================================================================== --- /branches/eam_branches/ipp-20121130/Ohana/src/relphot/include/relphot.h (revision 34836) +++ /branches/eam_branches/ipp-20121130/Ohana/src/relphot/include/relphot.h (revision 34837) @@ -63,4 +63,33 @@ typedef struct { + int Nfew; + int Ncode; + int Nsys; + int Nbad; + int Ncal; + int Nmos; + int Ngrid; + + // NOTE: the following arrays are (possibly) pre-allocated and carried down to each + // thread. The first 3 (list, dlist, wlist) are used in all relphot analyses; the + // others are only used on the final output steps. + + double *list; // list of measure.mag values for a given star + double *dlist; // mag errors for a star + double *wlist; // weights to use for mean mags + double *aplist; // ap mags for a star + double *kronlist; // kron mags for a star + double *dkronlist; // kron mag errors + + double *psfqf_list; // psfqf for all filters + double *psfqfperf_list; // psfqfperf for all filters + double *stargal_list; // stargal for all filters + + StatType stats; + StatType apstats; + StatType kronstats; +} SetMrelInfo; + +typedef struct { AverageTiny *average; // array of (minimal) average data MeasureTiny *measure; // array of (minimal) measure data @@ -113,4 +142,6 @@ char *BOUNDARY_TREE; + +int SET_MREL_VERSION; double MAG_LIM; @@ -239,4 +270,6 @@ Coords *getCoords PROTO((off_t meas, int cat)); off_t getImageEntry PROTO((off_t meas, int cat)); + +float getMcal_alt PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, float Xccd, float Yccd)); float getMcal PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog)); float getMflat PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog)); @@ -348,2 +381,10 @@ int load_tree (char *treefile); int BoundaryTreePrimaryCell (char *primaryCellName, double ra, double dec); + +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 setMrelAverageStack (off_t meas, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, char *primaryCell); +int setGlobalObjStats (Average *average, Measure *measure); + +void SetMrelInfoInit (SetMrelInfo *results, int allocLists); Index: /branches/eam_branches/ipp-20121130/Ohana/src/relphot/src/ImageOps.c =================================================================== --- /branches/eam_branches/ipp-20121130/Ohana/src/relphot/src/ImageOps.c (revision 34836) +++ /branches/eam_branches/ipp-20121130/Ohana/src/relphot/src/ImageOps.c (revision 34837) @@ -345,4 +345,27 @@ if (flat_id) { offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[cat].measureT[meas].Xccd, catalog[cat].measureT[meas].Yccd); + } + value -= offset; + + return (value); +} + +// returns image.Mcal - ff(x,y) +float getMcal_alt (off_t meas, int cat, FlatCorrectionTable *flatcorr, float Xccd, float Yccd) { + + off_t i; + float value, offset; + + i = MeasureToImage[cat][meas]; + if (i == -1) return (NAN); + + if (image[i].flags & IMAGE_BAD) return (NAN); + value = image[i].Mcal; + offset = 0.0; + + // to do this, I need to pass in the catalog and flatcorr pointers + int flat_id = image[i].photom_map_id; + if (flat_id) { + offset = FlatCorrectionOffset (flatcorr, flat_id, Xccd, Yccd); } value -= offset; Index: /branches/eam_branches/ipp-20121130/Ohana/src/relphot/src/StarOps.c =================================================================== --- /branches/eam_branches/ipp-20121130/Ohana/src/relphot/src/StarOps.c (revision 34836) +++ /branches/eam_branches/ipp-20121130/Ohana/src/relphot/src/StarOps.c (revision 34837) @@ -2,29 +2,4 @@ static int Nmax; - -typedef struct { - int Nfew; - int Ncode; - int Nsys; - int Nbad; - int Ncal; - int Nmos; - int Ngrid; - - // NOTE: the following arrays are (possibly) pre-allocated and carried down to each - // thread. The first 3 (list, dlist, wlist) are used in all relphot analyses; the - // others are only used on the final output steps. - - double *list; // list of measure.mag values for a given star - double *dlist; // mag errors for a star - double *wlist; // weights to use for mean mags - double *aplist; // ap mags for a star - double *kronlist; // kron mags for a star - double *dkronlist; // kron mag errors - - double *psfqf_list; // psfqf for all filters - double *psfqfperf_list; // psfqfperf for all filters - double *stargal_list; // stargal for all filters -} SetMrelInfo; enum {THREAD_RUN, THREAD_DONE}; @@ -154,5 +129,15 @@ for (i = 0; i < Ncatalog; i++) { // pass == -1 for anything other than the final pass - setMrel_catalog (catalog, i, -1, flatcorr, &results, Nsecfilt); + switch (SET_MREL_VERSION) { + case 0: + setMrel_catalog (catalog, i, -1, flatcorr, &results, Nsecfilt); + break; + case 1: + setMrel_catalog_alt (catalog, i, -1, flatcorr, &results, Nsecfilt); + break; + default: + fprintf (stderr, "invalid setMrel version %d (use 0 or 1)\n", SET_MREL_VERSION); + exit (2); + } SetMrelInfoAccum (&summary, &results); } @@ -182,5 +167,15 @@ for (i = 0; i < Ncatalog; i++) { - setMrel_catalog (catalog, i, pass, flatcorr, &results, Nsecfilt); // XXX add arguments as needed for options + switch (SET_MREL_VERSION) { + case 0: + setMrel_catalog (catalog, i, pass, flatcorr, &results, Nsecfilt); + break; + case 1: + setMrel_catalog_alt (catalog, i, pass, flatcorr, &results, Nsecfilt); + break; + default: + fprintf (stderr, "invalid setMrel version %d (use 0 or 1)\n", SET_MREL_VERSION); + exit (2); + } SetMrelInfoAccum (&summary, &results); } @@ -286,5 +281,15 @@ // pass == -1 for anything other than the final pass - setMrel_catalog (catalog, i, -1, flatcorr, &results, Nsecfilt); + switch (SET_MREL_VERSION) { + case 0: + setMrel_catalog (catalog, i, -1, flatcorr, &results, Nsecfilt); + break; + case 1: + setMrel_catalog_alt (catalog, i, -1, flatcorr, &results, Nsecfilt); + break; + default: + fprintf (stderr, "invalid setMrel version %d (use 0 or 1)\n", SET_MREL_VERSION); + exit (2); + } SetMrelInfoAccum (&threadinfo->summary, &results); } @@ -331,7 +336,4 @@ liststats_setmode (&apstats, STATMODE); liststats_setmode (&kronstats, STATMODE); - - // StatType psfstats; - // liststats_setmode (&psfstats, "MEAN"); double *list = results->list; Index: /branches/eam_branches/ipp-20121130/Ohana/src/relphot/src/args.c =================================================================== --- /branches/eam_branches/ipp-20121130/Ohana/src/relphot/src/args.c (revision 34836) +++ /branches/eam_branches/ipp-20121130/Ohana/src/relphot/src/args.c (revision 34837) @@ -64,4 +64,11 @@ } + SET_MREL_VERSION = 1; + if ((N = get_argument (argc, argv, "-set-mrel-version"))) { + remove_argument (N, &argc, argv); + SET_MREL_VERSION = atof(argv[N]); + remove_argument (N, &argc, argv); + } + VERBOSE = VERBOSE2 = FALSE; if ((N = get_argument (argc, argv, "-v"))) { Index: /branches/eam_branches/ipp-20121130/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- /branches/eam_branches/ipp-20121130/Ohana/src/relphot/src/setMrelCatalog.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/Ohana/src/relphot/src/setMrelCatalog.c (revision 34837) @@ -0,0 +1,588 @@ +# include "relphot.h" + +# define SKIP_THIS_MEAS(REASON) { \ + measureT[k].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ + if (measure) { \ + measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ + } \ + results->REASON ++; \ + continue; } + +# define SKIP_THIS_MEAS_STACK(REASON) { \ + measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ + continue; } + +int print_measure_set_alt (Average *average, SecFilt *secfilt, Measure *measure) { + + int Ns; + off_t k; + + int Nsecfilt = GetPhotcodeNsecfilt (); + + for (k = 0; k < average[0].Nmeasure; k++) { + fprintf (stderr, "meas: %08x\n", measure[k].dbFlags); + } + + for (Ns = 0; Ns < Nsecfilt; Ns++) { + fprintf (stderr, "secf: %08x\n", secfilt[Ns].flags); + } + return (TRUE); +} + +int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) { + + off_t j; + + liststats_setmode (&results->stats, STATMODE); + liststats_setmode (&results->apstats, STATMODE); + liststats_setmode (&results->kronstats, STATMODE); + + 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++) { + off_t m = catalog[Nc].average[j].measureOffset; + setMrelAverageExposure (m, Nc, pass, flatcorr, results, &catalog[Nc].average[j], &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt], &catalog[Nc].measure[m], &catalog[Nc].measureT[m], &catalog[Nc].found[Nsecfilt*j]); + + // only apply Stack operation on setMrelFinal in first pass + if (isSetMrelFinal && (pass == 0)) { + setMrelAverageStack (m, Nc, flatcorr, results, &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], &catalog[Nc].measure[m], primaryCell); + setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]); + } + } + if (primaryCell) free (primaryCell); + return (TRUE); +} + +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) { + + off_t k; + int N; + float Msys = 0, Mcal= 0, Mmos = 0, Mgrid = 0; + + double *list = results->list; + double *dlist = results->dlist; + double *wlist = results->wlist; + double *aplist = results->aplist; + double *kronlist = results->kronlist; + double *dkronlist = results->dkronlist; + + StatType *stats = &results->stats; + 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 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; + + N = 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 ++; + + 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); + } + + Msys = PhotSysTiny (&measureT[k], &averageT[0], &secfilt[0]); + if (isnan(Msys)) SKIP_THIS_MEAS(Nsys); + + int myUbercalDist = getUbercalDist(meas, cat); + minUbercalDist = MIN(minUbercalDist, myUbercalDist); + + 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; + } + + // gpc1 stack data + if ((measure[k].photcode >= 11000) && (measure[k].photcode <= 11400)) { + continue; + } + + // 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 = 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 (&measure[k]); + aplist[N] = Map - Mcal - Mmos - Mgrid; + + float Mkron = PhotKron (&measure[k]); + kronlist[N] = Mkron - Mcal - Mmos - Mgrid; + dkronlist[N] = measure[k].dMkron; + } + + // 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(measureT[k].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 (measureT[k].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; + + // 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", averageT[0].R, averageT[0].D, measureT[0].imageID, N, 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 < 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); + + secfilt[Nsec].M = stats->mean; + secfilt[Nsec].dM = stats->error; + secfilt[Nsec].Xm = (stats->Nmeas > 1) ? 100.0*log10(stats->chisq + 1e-4) : NAN_S_SHORT; + + // 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 = 1000.0*stats->sigma; // Mstdev is in millimags (not enough space for more precision) + secfilt[Nsec].Ncode = Ncode; + secfilt[Nsec].Nused = stats->Nmeas; + + secfilt[Nsec].M_80 = 1000 * stats->Upper80; + secfilt[Nsec].M_20 = 1000 * stats->Lower20; + + // NOTE : use the modified weight for apmags as well as psf mags + liststats (aplist, dlist, wlist, N, apstats); + secfilt[Nsec].Map = apstats->mean; + + liststats (kronlist, dkronlist, wlist, N, kronstats); + secfilt[Nsec].Mkron = kronstats->mean; + secfilt[Nsec].dMkron = kronstats->error; + + // NOTE: for 2MASS measurements, Next should be 1, as should N + if ((Next > 0) && (Next > 0.5*N)) { + 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 +int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, char *primaryCell) { + + off_t k, ID; + + float Msys = 0, Mcal= 0, Mmos = 0, Mgrid = 0; + + // set the name of the primary skycell (this is used in a strcmp to match the skycells in stack detections) + BoundaryTreePrimaryCell(primaryCell, average[0].R, average[0].D); + + 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 */ + + int haveStack = FALSE; + + // need to find the measurement closest to the center of its skycell, as well as the + // closest for the subset of primary projection cells + + float stackCenterOffsetMin = 1e9; + off_t stackCenterMeasureMin = -1; + + float stackPrimaryOffsetMin = 1e9; + off_t stackPrimaryMeasureMin = -1; + + int isBad = FALSE; + int isSuspect = FALSE; + + off_t meas = measureOffset; + for (k = 0; k < average[0].Nmeasure; k++, meas++) { + + // skip measurements that do not match the current photcode + PhotCode *code = GetPhotcodebyCode (measure[k].photcode); + if (!code) continue; + if (code->equiv != thisCode) { continue; } + + // clear this bit for all measurements + measure[k].dbFlags &= ~ID_MEAS_STACK_PRIMARY; + + // only examine gpc1 stack data + // XXX this is absurdly hardwired (along with several photcode tests) + if (measure[k].photcode < 11000) continue; + if (measure[k].photcode > 11400) continue; + + haveStack = TRUE; + haveStackObject = TRUE; + + if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad); + + if (getImageEntry (meas, cat) < 0) continue; + // measurements without an image are either external reference photometry or + // data for which the associated image has not been loaded (probably because of + // overlaps). we only want measurements associated with stack images in this loop + + Mcal = getMcal_alt (meas, cat, flatcorr, measure[k].Xccd, measure[k].Yccd); + if (isnan(Mcal)) SKIP_THIS_MEAS_STACK(Ncal); + Mmos = getMmos (meas, cat); + if (isnan(Mmos)) SKIP_THIS_MEAS_STACK(Nmos); + Mgrid = getMgrid (meas, cat); + if (isnan(Mgrid)) SKIP_THIS_MEAS_STACK(Ngrid); + + Msys = PhotSys (&measure[k], &average[0], &secfilt[0]); + if (isnan(Msys)) SKIP_THIS_MEAS_STACK(Nsys); + + 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 (meas, cat, primaryCell)) { + float stackPrimaryOffset = getCenterOffset (meas, cat, &measure[k], &stackImageID); + if (stackPrimaryOffset < stackPrimaryOffsetMin) { + stackPrimaryOffsetMin = stackPrimaryOffset; + // stackPrimaryIDmin = stackImageID; + stackPrimaryMeasureMin = k; + } + } + + // get the center distance for the generic case: + float stackCenterOffset = getCenterOffset (meas, cat, &measure[k], &stackImageID); + if (stackCenterOffset < stackCenterOffsetMin) { + stackCenterOffsetMin = stackCenterOffset; + // stackCenterIDmin = stackImageID; + stackCenterMeasureMin = k; + } + + // data quality assessment + isBad |= (measure[k].photFlags & code->photomBadMask); + isBad |= (measure[k].psfQF < 0.85); + isBad |= isnan(measure[k].psfQF); + isBad |= measure[k].dM > 0.2; // S/N < 5.0 + + isSuspect |= (measure[k].photFlags & code->photomPoorMask); + isSuspect |= (measure[k].psfQFperf < 0.85); + } + + if (!haveStack) continue; + + if (!isSuspect && !isBad) { + NstackGood ++; + } + if (isSuspect && !isBad) { + NstackSuspect ++; + } + + // measurements which are bad will not have a valid stack entry and are skipped + k = (stackPrimaryMeasureMin >= 0) ? stackPrimaryMeasureMin : stackCenterMeasureMin; + if (k < 0) continue; + + + // we are now populating stackDetectID not stack Image ID in secfilt + // ID = (stackPrimaryMeasureMin >= 0) ? stackPrimaryIDmin : stackCenterIDmin; + ID = measure[k].extID; // for the stack, this is the stackDetectID + + // get the zero point for the selected image + float zp = PhotZeroPoint (&measure[k], &average[0], &secfilt[0]) - (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?) + secfilt[Nsec].FluxPSF = zpFactor * measure[k].FluxPSF; + secfilt[Nsec].dFluxPSF = zpFactor * measure[k].dFluxPSF; + secfilt[Nsec].FluxKron = zpFactor * measure[k].FluxKron; + secfilt[Nsec].dFluxKron = zpFactor * measure[k].dFluxKron; + + secfilt[Nsec].stackDetectID = ID; + + measure[k].dbFlags |= ID_MEAS_STACK_PHOT_SRC; + if (stackPrimaryMeasureMin >= 0) { + measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY; + } + + secfilt[Nsec].flags |= ID_SECF_HAS_STACK; + } // Nsecfilt loop + + if (NstackGood >= 2) { + average[0].flags |= ID_OBJ_GOOD_STACK; + } else if (NstackGood + NstackSuspect >= 2) { + average[0].flags |= ID_OBJ_SUSPECT_STACK; + } else if (haveStackObject) { + average[0].flags |= ID_OBJ_BAD_STACK; + } + return (TRUE); +} + +int setGlobalObjStats (Average *average, Measure *measure) { + + off_t k; + + float psfQFmax = 0.0; + float psfQFperfmax = 0.0; + float stargalmax = 0.0; + + for (k = 0; k < average[0].Nmeasure; k++) { + + if (isfinite(measure[k].psfQF)) { + psfQFmax = MAX (psfQFmax, measure[k].psfQF); + } + if (isfinite(measure[k].psfQFperf)) { + psfQFperfmax = MAX (psfQFperfmax, measure[k].psfQFperf); + } + + // only count psfQF > 0.85 here + if (isfinite(measure[k].extNsigma) && (measure[k].psfQF > 0.85)) { + stargalmax = MAX (stargalmax, measure[k].extNsigma); + } + } + average[0].psfQF = psfQFmax; + average[0].psfQFperf = psfQFperfmax; + average[0].stargal = stargalmax; + + return (TRUE); +}