Index: trunk/Ohana/src/relphot/src/StarOps.c =================================================================== --- trunk/Ohana/src/relphot/src/StarOps.c (revision 34405) +++ trunk/Ohana/src/relphot/src/StarOps.c (revision 34429) @@ -11,10 +11,19 @@ int Nmos; int Ngrid; - double *list; - double *dlist; - double *wlist; - double *aplist; - double *kronlist; - double *dkronlist; + + // 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; @@ -168,4 +177,8 @@ ALLOCATE (results.dkronlist, double, Nmax); + ALLOCATE (results.psfqf_list, double, Nmax); + ALLOCATE (results.psfqfperf_list, double, Nmax); + ALLOCATE (results.stargal_list, double, Nmax); + for (i = 0; i < Ncatalog; i++) { setMrel_catalog (catalog, i, pass, flatcorr, &results, Nsecfilt); // XXX add arguments as needed for options @@ -178,4 +191,8 @@ free (results.kronlist); free (results.dkronlist); + + free (results.psfqf_list); + free (results.psfqfperf_list); + free (results.stargal_list); return (TRUE); } @@ -315,4 +332,7 @@ liststats_setmode (&kronstats, STATMODE); + // StatType psfstats; + // liststats_setmode (&psfstats, "MEAN"); + double *list = results->list; double *dlist = results->dlist; @@ -322,4 +342,8 @@ 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 @@ -327,5 +351,7 @@ char *primaryCell = NULL; - ALLOCATE (primaryCell, char, DVO_MAX_PATH); + if (isSetMrelFinal) { + ALLOCATE (primaryCell, char, DVO_MAX_PATH); + } for (j = 0; j < catalog[Nc].Naverage; j++) { @@ -338,5 +364,7 @@ } - BoundaryTreePrimaryCell(primaryCell, catalog[Nc].average[j].R, catalog[Nc].average[j].D); + if (isSetMrelFinal) { + BoundaryTreePrimaryCell(primaryCell, catalog[Nc].average[j].R, catalog[Nc].average[j].D); + } int GoodPS1 = FALSE; @@ -346,4 +374,6 @@ int NextPS1 = 0; int NpsfPS1 = 0; + + int Nmeas = 0; // number of meas used for psfqf, psfqfperf, and stargal int Ns; @@ -365,4 +395,5 @@ 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 @@ -420,4 +451,5 @@ if (isSetMrelFinal) { + // 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; @@ -426,4 +458,9 @@ 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) + psfqf_list[Nmeas] = catalog[Nc].measure[m].psfQF; + psfqfperf_list[Nmeas] = catalog[Nc].measure[m].psfQF; + stargal_list[Nmeas] = catalog[Nc].measure[m].extNsigma; // special options for PS1 data @@ -431,5 +468,10 @@ // count the extended detections if (!isnan(catalog[Nc].measure[m].Map)) { - if (catalog[Nc].measure[m].M - catalog[Nc].measure[m].Map > 0.5) { + 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 ++; @@ -438,4 +480,5 @@ } } + havePS1 = TRUE; } // gpc1 stack data @@ -543,4 +586,5 @@ } N++; + Nmeas++; } @@ -658,13 +702,13 @@ case 0: catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_0; - GoodPS1 = TRUE; + if (havePS1) GoodPS1 = TRUE; break; case 1: catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_1; - GoodPS1 = TRUE; + if (havePS1) GoodPS1 = TRUE; break; case 2: catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_2; - GoodPS1 = TRUE; + if (havePS1) GoodPS1 = TRUE; break; case 3: @@ -678,4 +722,10 @@ 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; @@ -685,8 +735,9 @@ if (isSetMrelFinal) { - 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 (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)) { @@ -701,4 +752,28 @@ if (Good2MASS) { catalog[Nc].average[j].flags |= ID_OBJ_GOOD_ALT; + } + + // 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; } } Index: trunk/Ohana/src/relphot/src/initialize.c =================================================================== --- trunk/Ohana/src/relphot/src/initialize.c (revision 34405) +++ trunk/Ohana/src/relphot/src/initialize.c (revision 34429) @@ -77,5 +77,21 @@ args_client (argc, argv); - if (MODE == MODE_UPDATE_OBJECTS) return; + if (MODE == MODE_UPDATE_OBJECTS) { + char tmpline1[256]; + int Ns; + Nphotcodes = GetPhotcodeNsecfilt (); + ALLOCATE (photcodes, PhotCode *, Nphotcodes); + ALLOCATE (PhotcodeList, char, 256); + for (Ns = 0; Ns < Nphotcodes; Ns++) { + photcodes[Ns] = GetPhotcodebyNsec (Ns); + if (Ns > 0) { + snprintf (tmpline1, 256, "%s,%s", PhotcodeList, photcodes[Ns][0].name); + } else { + snprintf (tmpline1, 256, "%s", photcodes[Ns][0].name); + } + strcpy (PhotcodeList, tmpline1); + } + return; + } // load the list of photcodes into the globals (photcodes, Nphotcodes) Index: trunk/Ohana/src/relphot/src/liststats.c =================================================================== --- trunk/Ohana/src/relphot/src/liststats.c (revision 34405) +++ trunk/Ohana/src/relphot/src/liststats.c (revision 34429) @@ -38,5 +38,9 @@ dsortthree (value, dvalue, weight, N); } else { - dsortpair (value, dvalue, N); + if (dvalue) { + dsortpair (value, dvalue, N); + } else { + dsort (value, N); + } } Index: trunk/Ohana/src/relphot/src/relphot_objects.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_objects.c (revision 34405) +++ trunk/Ohana/src/relphot/src/relphot_objects.c (revision 34429) @@ -61,29 +61,34 @@ Nsecfilt = catalog.Nsecfilt; DVOSecfiltFlags secfiltBits = - ID_SECF_STAR_FEW | - ID_SECF_STAR_POOR | - ID_PHOTOM_PASS_0 | - ID_PHOTOM_PASS_1 | - ID_PHOTOM_PASS_2 | - ID_PHOTOM_PASS_3 | - ID_PHOTOM_PASS_4 | - ID_SECF_USE_SYNTH | + ID_SECF_STAR_FEW | + ID_SECF_STAR_POOR | + ID_PHOTOM_PASS_0 | + ID_PHOTOM_PASS_1 | + ID_PHOTOM_PASS_2 | + ID_PHOTOM_PASS_3 | + ID_PHOTOM_PASS_4 | + ID_SECF_USE_SYNTH | ID_SECF_USE_UBERCAL | - ID_SECF_OBJ_EXT; + ID_SECF_HAS_STACK | + ID_SECF_HAS_PS1 | + ID_SECF_OBJ_EXT ; + + DVOAverageFlags astromBits = + ID_STAR_FIT_AVE | + ID_STAR_FIT_PM | + ID_STAR_FIT_PAR | + ID_STAR_USE_AVE | + ID_STAR_USE_PM | + ID_STAR_USE_PAR | + ID_STAR_NO_ASTROM ; for (j = 0; j < catalog.Naverage; j++) { - catalog.average[j].flags = 0; + catalog.average[j].flags &= astromBits; // 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) for (k = 0; k < Nsecfilt; k++) { - catalog.secfilt[j*Nsecfilt + k].M = NAN; - catalog.secfilt[j*Nsecfilt + k].Map = NAN; - catalog.secfilt[j*Nsecfilt + k].dM = NAN; - catalog.secfilt[j*Nsecfilt + k].Mstdev = NAN_S_SHORT; - catalog.secfilt[j*Nsecfilt + k].Xm = NAN; - catalog.secfilt[j*Nsecfilt + k].M_20 = NAN_S_SHORT; - catalog.secfilt[j*Nsecfilt + k].M_80 = NAN_S_SHORT; - catalog.secfilt[j*Nsecfilt + k].Ncode = 0; - catalog.secfilt[j*Nsecfilt + k].Nused = 0; - catalog.secfilt[j*Nsecfilt + k].flags &= ~secfiltBits; - catalog.secfilt[j*Nsecfilt + k].ubercalDist = 1000; + dvo_secfilt_init (&catalog.secfilt[j*Nsecfilt + k]); + catalog.secfilt[j*Nsecfilt + k].flags &= ~secfiltBits; // XXX: only reset photom bits: } } Index: trunk/Ohana/src/relphot/src/setMrelFinal.c =================================================================== --- trunk/Ohana/src/relphot/src/setMrelFinal.c (revision 34405) +++ trunk/Ohana/src/relphot/src/setMrelFinal.c (revision 34429) @@ -94,4 +94,10 @@ } + for (i = 0; i < catalog[0].Naverage; i++) { + catalog[0].average[i].psfQF = NAN; // force recalculation below + catalog[0].average[i].psfQFperf = NAN; // force recalculation below + catalog[0].average[i].stargal = NAN; // force recalculation below + } + // XXX make this optional? (do not clean for -averages?) if (!simpleAverage) clean_measures (catalog, 1, TRUE, flatcorr); /* mark outliers ID_MEAS_POOR_PHOTOM */