Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/include/relphot.h =================================================================== --- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/include/relphot.h (revision 33822) +++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/include/relphot.h (revision 33823) @@ -277,5 +277,5 @@ int setMmos PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr)); int setMrel PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); -void setMrelFinal PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); +void setMrelFinal PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage)); int setMrelOutput PROTO((Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr)); int setMave PROTO((Catalog *catalog, int Ncatalog)); Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/ImageOps.c =================================================================== --- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/ImageOps.c (revision 33822) +++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/ImageOps.c (revision 33823) @@ -41,7 +41,7 @@ // relationships between the measure,catalog set and the images: -static IDX_T **MeasureToImage; // image index from measure,catalog : MeasureToImage[cat][meas] = ImageIndex -static IDX_T **ImageToCatalog; // catalog for given measure on image : ImageCatalog[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex]) -static IDX_T **ImageToMeasure; // measure for given measure on image : ImageMeasure[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex]) +static IDX_T **MeasureToImage = NULL; // image index from measure,catalog : MeasureToImage[cat][meas] = ImageIndex +static IDX_T **ImageToCatalog = NULL; // catalog for given measure on image : ImageCatalog[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex]) +static IDX_T **ImageToMeasure = NULL; // measure for given measure on image : ImageMeasure[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex]) // MeasureToImage was 'bin' @@ -322,4 +322,6 @@ off_t i; + if (!MeasureToImage) return -1; + i = MeasureToImage[cat][meas]; return (i); @@ -354,4 +356,6 @@ off_t i; short distance; + + if (!MeasureToImage) return -1; i = MeasureToImage[cat][meas]; Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/StarOps.c =================================================================== --- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/StarOps.c (revision 33822) +++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/StarOps.c (revision 33823) @@ -5,4 +5,5 @@ typedef struct { int Nfew; + int Ncode; int Nsys; int Nbad; @@ -12,4 +13,6 @@ double *list; double *dlist; + double *aplist; + double *daplist; } SetMrelInfo; @@ -25,7 +28,9 @@ } ThreadInfo; +// ** internal functions: void *setMrel_worker (void *data); int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr); -int setMrel_catalog (Catalog *catalog, int Nc, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt); +int setMrel_catalog (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt); +int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure); // we want to allocate the stats list,dlist arrays only once (or once per thread). @@ -65,4 +70,5 @@ void SetMrelInfoInit (SetMrelInfo *results, int allocLists) { results->Nfew = 0; + results->Ncode = 0; results->Nsys = 0; results->Nbad = 0; @@ -83,4 +89,5 @@ void SetMrelInfoAccum (SetMrelInfo *summary, SetMrelInfo *results) { summary->Nfew += results->Nfew ; + summary->Ncode += results->Ncode ; summary->Nsys += results->Nsys ; summary->Nbad += results->Nbad ; @@ -110,4 +117,11 @@ } +// setMrel and setMrelOutput are extremely similar, but have slightly different implications: +// * setMrel uses the internal Tiny structures only +// * setMrelOutput skips stars for which there are too few good measurements +// * setMrelOutput is meant to be called repeatedly, relaxing the criteria for 'good' on each pass +// * setMrelOutput updates 2MASS average flags +// * setMrelOutput updates average EXT flags (PS1 and 2MASS) + int setMrel (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { @@ -126,8 +140,9 @@ for (i = 0; i < Ncatalog; i++) { - setMrel_catalog (catalog, i, flatcorr, &results, Nsecfilt); + // pass == -1 for anything other than the final pass + setMrel_catalog (catalog, i, -1, flatcorr, &results, Nsecfilt); SetMrelInfoAccum (&summary, &results); } - fprintf (stderr, "%d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys); + fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys); SetMrelInfoFree (&results); @@ -135,23 +150,179 @@ return (TRUE); } + +int setMrelOutput (Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr) { + + int i; + + int Nsecfilt = GetPhotcodeNsecfilt (); + + SetMrelInfo summary, results; + SetMrelInfoInit (&summary, FALSE); + SetMrelInfoInit (&results, TRUE); // allocates results->list,dlist + ALLOCATE (results.aplist, double, Nmax); + ALLOCATE (results.daplist, double, Nmax); + + for (i = 0; i < Ncatalog; i++) { + setMrel_catalog (catalog, i, pass, flatcorr, &results, Nsecfilt); // XXX add arguments as needed for options + SetMrelInfoAccum (&summary, &results); + } + fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys); + + SetMrelInfoFree (&results); + free (results.aplist); + free (results.daplist); + return (TRUE); +} + +int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { + + int i; + + SetMrelInfo summary; + SetMrelInfoInit (&summary, FALSE); + + pthread_attr_t attr; + pthread_attr_init (&attr); + pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED); + + pthread_t *threads; + ALLOCATE (threads, pthread_t, NTHREADS); + + ThreadInfo *threadinfo; + ALLOCATE (threadinfo, ThreadInfo, NTHREADS); + + // launch N worker threads + for (i = 0; i < NTHREADS; i++) { + threadinfo[i].entry = i; + threadinfo[i].state = THREAD_RUN; + threadinfo[i].catalog = catalog; + threadinfo[i].Ncatalog = Ncatalog; + threadinfo[i].flatcorr = flatcorr; + SetMrelInfoInit (&threadinfo[i].summary, FALSE); + pthread_create (&threads[i], NULL, setMrel_worker, &threadinfo[i]); + } + pthread_attr_destroy (&attr); + + // wait until all threads have finished + while (1) { + int allDone = TRUE; + for (i = 0; i < NTHREADS; i++) { + if (threadinfo[i].state == THREAD_RUN) allDone = FALSE; + } + if (allDone) { + break; + } + usleep (500000); + } + + // all threads are done, free the threads array and grab the info + free (threads); + + // report stats & summary from the threads + for (i = 0; i < NTHREADS; i++) { + fprintf (stderr, "setMrel thread %d : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", + i, + threadinfo[i].summary.Ncode, + threadinfo[i].summary.Nfew, + threadinfo[i].summary.Nbad, + threadinfo[i].summary.Ncal, + threadinfo[i].summary.Nmos, + threadinfo[i].summary.Ngrid, + threadinfo[i].summary.Nsys); + SetMrelInfoAccum (&summary, &threadinfo[i].summary); + } + fprintf (stderr, "total : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys); + free (threadinfo); + + return TRUE; +} + +void *setMrel_worker (void *data) { + + ThreadInfo *threadinfo = data; + + int Nsecfilt = GetPhotcodeNsecfilt (); + + SetMrelInfo results; + SetMrelInfoInit (&results, TRUE); // allocate list, dlist arrays here + + while (1) { + + off_t i = getNextCatalogForThread(threadinfo->Ncatalog); + if (i == -1) { + threadinfo->state = THREAD_DONE; + return NULL; + } + + Catalog *catalog = threadinfo->catalog; + FlatCorrectionTable *flatcorr = threadinfo->flatcorr; + + // pass == -1 for anything other than the final pass + setMrel_catalog (catalog, i, -1, flatcorr, &results, Nsecfilt); + SetMrelInfoAccum (&threadinfo->summary, &results); + } + + SetMrelInfoFree (&results); + return NULL; +} + +# define SKIP_THIS_MEAS(REASON) { \ + catalog[Nc].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ + if (catalog[Nc].measure) { \ + catalog[Nc].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ + } \ + results->REASON ++; \ + 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, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) { +int setMrel_catalog (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) { off_t j, k, m; int N; float Msys, Mcal, Mmos, Mgrid; - StatType stats; - - double *list = results->list; - double *dlist = results->dlist; + StatType stats, apstats; + + double *list = results->list; + double *dlist = results->dlist; + double *aplist = results->aplist; + double *daplist = results->daplist; SetMrelInfoInit (results, FALSE); // do not allocate list,dlist arrays + + int isSetMrelFinal = (pass >= 0); for (j = 0; j < catalog[Nc].Naverage; j++) { // XXX accumulate all secfilt values in a single pass? - int minUbercalDist = 1000; + // option for a test print + if (FALSE && (catalog[Nc].average[j].objID == 0x46a4) && (catalog[Nc].average[j].catID == 0xf40e)) { + fprintf (stderr, "test obj\n"); + print_measure_set (&catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], catalog[Nc].measure); + } + + int GoodPS1 = FALSE; + int Good2MASS = FALSE; + int Galaxy2MASS = FALSE; + + int NextPS1 = 0; + int NpsfPS1 = 0; int Ns; @@ -163,7 +334,17 @@ /* calculate the average mag in this SEC photcode for a single star */ + /* star/photcodes already calibrated */ + if ( isSetMrelFinal && catalog[Nc].found[Nsecfilt*j+Nsec]) continue; + // skip bad stars - if (catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; - + if (!isSetMrelFinal && (catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD)) continue; + + int Ncode = 0; + int Next = 0; + int haveSynth = FALSE; + int haveUbercal = FALSE; + + int minUbercalDist = 1000; + N = 0; m = catalog[Nc].averageT[j].measureOffset; @@ -174,6 +355,7 @@ if (!code) continue; if (code->equiv != thisCode) { continue; } - - if (catalog[Nc].measureT[m].dbFlags & MEAS_BAD) { results->Nbad ++; continue; } + Ncode ++; + + if (catalog[Nc].measureT[m].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); if (getImageEntry (m, Nc) < 0) { @@ -185,17 +367,59 @@ } else { Mcal = getMcal (m, Nc, flatcorr, catalog); - if (isnan(Mcal)) { results->Ncal ++; continue; } + if (isnan(Mcal)) SKIP_THIS_MEAS(Ncal); Mmos = getMmos (m, Nc); - if (isnan(Mmos)) { results->Nmos ++; continue; } + if (isnan(Mmos)) SKIP_THIS_MEAS(Nmos); Mgrid = getMgrid (m, Nc); - if (isnan(Mgrid)) { results->Ngrid++; continue; } + if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid); } Msys = PhotSysTiny (&catalog[Nc].measureT[m], &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt]); - if (isnan(Msys)) { results->Nsys++; continue; } + if (isnan(Msys)) SKIP_THIS_MEAS(Nsys); + list[N] = Msys - Mcal - Mmos - Mgrid; int myUbercalDist = getUbercalDist(m,Nc); minUbercalDist = MIN(minUbercalDist, myUbercalDist); + + if (isSetMrelFinal) { + float Map = PhotAper (&catalog[Nc].measure[m]); + aplist[N] = Map - Mcal - Mmos - Mgrid; + + // count the extended detections + if ((catalog[Nc].measure[m].photcode >= 10000) && (catalog[Nc].measure[m].photcode <= 10500)) { + if (!isnan(catalog[Nc].measure[m].Map)) { + if (catalog[Nc].measure[m].M - catalog[Nc].measure[m].Map > 0.5) { + Next ++; + NextPS1 ++; + } else { + NpsfPS1 ++; + } + } + } + // 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); + continue; + } + } + } + // ignore SYNTH photocdes until PASS == 4 (where we also accept saturated stars) + if ((catalog[Nc].measure[m].photcode >= 3001) && (catalog[Nc].measure[m].photcode <= 3005)) { + if (pass < 4) { + SKIP_THIS_MEAS(Nbad); + continue; + } + haveSynth = TRUE; + } + } // dlist gives the error, which is used as the weight in WT_MEAN. @@ -225,9 +449,16 @@ } - // when performing the grid analysis, STAR_TOOFEW will be set to 1; - if (N <= STAR_TOOFEW) { /* too few measurements */ + 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", 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_STAR_FEW; - results->Nfew ++; + if (Ncode == 0) { + results->Ncode ++; + } else { + results->Nfew ++; + } + continue; } else { catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_FEW; @@ -235,105 +466,80 @@ liststats (list, dlist, N, &stats); - + catalog[Nc].secfilt[Nsecfilt*j+Nsec].M = stats.mean; - catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM = stats.sigma; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM = stats.error; catalog[Nc].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq + 1e-4) : NAN_S_SHORT; - - catalog[Nc].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist; - } - } - - // values which need to be passed out: Nfew, Nbad, Ncal, Nmos, Ngrid, Nsys + + // 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) { + liststats (aplist, daplist, N, &apstats); + catalog[Nc].found[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: 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; + GoodPS1 = TRUE; + break; + case 1: + catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_1; + GoodPS1 = TRUE; + break; + case 2: + catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_2; + 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 (haveUbercal) { + catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_UBERCAL; + } + } + } + + 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 (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; + } + } + } return (TRUE); -} - -int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { - - int i; - - SetMrelInfo summary; - SetMrelInfoInit (&summary, FALSE); - - pthread_attr_t attr; - pthread_attr_init (&attr); - pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED); - - pthread_t *threads; - ALLOCATE (threads, pthread_t, NTHREADS); - - ThreadInfo *threadinfo; - ALLOCATE (threadinfo, ThreadInfo, NTHREADS); - - // launch N worker threads - for (i = 0; i < NTHREADS; i++) { - threadinfo[i].entry = i; - threadinfo[i].state = THREAD_RUN; - threadinfo[i].catalog = catalog; - threadinfo[i].Ncatalog = Ncatalog; - threadinfo[i].flatcorr = flatcorr; - SetMrelInfoInit (&threadinfo[i].summary, FALSE); - pthread_create (&threads[i], NULL, setMrel_worker, &threadinfo[i]); - } - pthread_attr_destroy (&attr); - - // wait until all threads have finished - while (1) { - int allDone = TRUE; - for (i = 0; i < NTHREADS; i++) { - if (threadinfo[i].state == THREAD_RUN) allDone = FALSE; - } - if (allDone) { - break; - } - usleep (500000); - } - - // all threads are done, free the threads array and grab the info - free (threads); - - // report stats & summary from the threads - for (i = 0; i < NTHREADS; i++) { - fprintf (stderr, "setMrel thread %d : %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", - i, - threadinfo[i].summary.Nfew, - threadinfo[i].summary.Nbad, - threadinfo[i].summary.Ncal, - threadinfo[i].summary.Nmos, - threadinfo[i].summary.Ngrid, - threadinfo[i].summary.Nsys); - SetMrelInfoAccum (&summary, &threadinfo[i].summary); - } - fprintf (stderr, "total : %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys); - free (threadinfo); - - return TRUE; -} - -void *setMrel_worker (void *data) { - - ThreadInfo *threadinfo = data; - - int Nsecfilt = GetPhotcodeNsecfilt (); - - SetMrelInfo results; - SetMrelInfoInit (&results, TRUE); // allocate list, dlist arrays here - - while (1) { - - off_t i = getNextCatalogForThread(threadinfo->Ncatalog); - if (i == -1) { - threadinfo->state = THREAD_DONE; - return NULL; - } - - Catalog *catalog = threadinfo->catalog; - FlatCorrectionTable *flatcorr = threadinfo->flatcorr; - - setMrel_catalog (catalog, i, flatcorr, &results, Nsecfilt); - SetMrelInfoAccum (&threadinfo->summary, &results); - } - - SetMrelInfoFree (&results); - return NULL; } @@ -357,233 +563,5 @@ } -# define MARK_SKIP_MEAS \ - catalog[i].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ - catalog[i].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM; - -// setMrel and setMrelOutput are extremely similar, but have slightly different implications: -// * setMrel uses the internal Tiny structures only -// * setMrelOutput skips stars for which there are too few good measurements -// * setMrelOutput is meant to be called repeatedly, relaxing the criteria for 'good' on each pass - -// setMave is also similar to the above. but: -// * setMave is called by relphot -update-objects -// * setMave excludes all detections with (PSF_QF < 0.85), setMrelOutput allows these for PASS > 2 -// * setMave updates 2MASS average flags -// * setMave updates average EXT flags (PS1 and 2MASS) -int setMrelOutput (Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr) { - - int i, N; - off_t j, k, m, Nmax; - float Msys, Mcal, Mmos, Mgrid, Map; - double *list, *dlist, *aplist, *daplist; - StatType stats, apstats; - int Nsec; - - int Nsecfilt = GetPhotcodeNsecfilt (); - - /* Nmeasure is now different, need to reallocate */ - Nmax = 0; - for (i = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { - Nmax = MAX (Nmax, catalog[i].averageT[j].Nmeasure); - } - } - ALLOCATE (list, double, MAX (1, Nmax)); - ALLOCATE (dlist, double, MAX (1, Nmax)); - - ALLOCATE (aplist, double, MAX (1, Nmax)); - ALLOCATE (daplist, double, MAX (1, Nmax)); - - for (i = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { - - if (FALSE && (catalog[i].average[j].objID == 0x46a4) && (catalog[i].average[j].catID == 0xf40e)) { - fprintf (stderr, "test obj\n"); - print_measure_set (&catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt], catalog[i].measure); - } - - int minUbercalDist = 1000; - - int Ns; - for (Ns = 0; Ns < Nphotcodes; Ns++) { - - int thisCode = photcodes[Ns][0].code; - Nsec = GetPhotcodeNsec(thisCode); - - /* star/photcodes already calibrated */ - if (catalog[i].found[Nsecfilt*j+Nsec]) continue; - - int Ncode = 0; - int Next = 0; - int haveSynth = FALSE; - int haveUbercal = FALSE; - int isRefPhot = FALSE; - - N = 0; - m = catalog[i].averageT[j].measureOffset; - for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { - - // skip measurements that do not match the current photcode - PhotCode *code = GetPhotcodebyCode (catalog[i].measureT[m].photcode); - if (!code) continue; - if (code->equiv != thisCode) { continue; } - Ncode ++; - - if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; - - if (getImageEntry (m, i) < 0) { - XXXX : likely problem -- compare with setMrel_catalog; - // these detetions have no image (eg, ref values such as 2MASS) - Mcal = Mmos = Mgrid = 0; - isRefPhot = TRUE; - } else { - Mcal = getMcal (m, i, flatcorr, catalog); - if (isnan(Mcal)) { - MARK_SKIP_MEAS; - continue; - } - Mmos = getMmos (m, i); - if (isnan(Mmos)) { - MARK_SKIP_MEAS; - continue; - } - Mgrid = getMgrid (m, i); - if (isnan(Mgrid)) { - MARK_SKIP_MEAS; - continue; - } - } - - Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]); - if (isnan(Msys)) { - if (!isRefPhot) { - MARK_SKIP_MEAS; - } - continue; - } - list[N] = Msys - Mcal - Mmos - Mgrid; - - Map = PhotAper (&catalog[i].measure[m]); - aplist[N] = Map - Mcal - Mmos - Mgrid; - - int myUbercalDist = getUbercalDist(m,i); - minUbercalDist = MIN(minUbercalDist, myUbercalDist); - - // count the extended detections - if ((catalog[i].measure[m].photcode >= 10000) && (catalog[i].measure[m].photcode <= 10500)) { - if (!isnan(catalog[i].measure[m].Map)) { - if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) { - Next ++; - } - } - } - // count extended detections for 2MASS (XXX NOTE hardwired photcodes 2011, 2012, 2013) - if ((catalog[i].measure[m].photcode >= 2011) && (catalog[i].measure[m].photcode <= 2013)) { - if (catalog[i].measure[m].photFlags & 0x00c00000) { - Next ++; - } - if ((pass == 0) && !(catalog[i].measure[m].photFlags & 0x00000007)) { - // detections without one of these bits should only be used in PASS_1 - MARK_SKIP_MEAS; - continue; - } - } - // ignore SYNTH photocdes until PASS == 4 (where we also accept saturated stars) - if ((catalog[i].measure[m].photcode >= 3001) && (catalog[i].measure[m].photcode <= 3005)) { - if (pass < 4) { - MARK_SKIP_MEAS; - continue; - } - haveSynth = TRUE; - } - - // dlist gives the error, which is used as the weight in WT_MEAN. - // we can modify the resulting 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[i].measureT[m].dM, code->photomErrSys), MIN_ERROR); - - // up-weight the ubercal values (or convergence can take a long time...) (XXX make this optional?) - if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) { - dlist[N] = MAX (0.1*catalog[i].measureT[m].dM, MIN_ERROR); - haveUbercal = TRUE; - } - - // 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) { - // tiny error -> large weight - // dlist[N] = MAX (0.01*catalog[i].measureT[m].dM, MIN_ERROR); - dlist[N] = 0.0001; - } - } - N++; - } - if (N < 1) continue; - - for (m = 0; m < N; m++) { - daplist[m] = dlist[m]; - } - - // XXX force WT_MEAN or MEAN here? - liststats (list, dlist, N, &stats); - liststats (aplist, daplist, N, &apstats); - catalog[i].found[Nsecfilt*j+Nsec] = TRUE; - switch (pass) { - case 0: - catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_0; - break; - case 1: - catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_1; - break; - case 2: - catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_2; - break; - case 3: - catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_3; - break; - case 4: - catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_4; - break; - } - - if (haveSynth) { - catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_SYNTH; - } - if (haveUbercal) { - catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_UBERCAL; - } - - /* use sigma or error in dM for output? */ - catalog[i].secfilt[Nsecfilt*j+Nsec].M = stats.mean; - catalog[i].secfilt[Nsecfilt*j+Nsec].Map = apstats.mean; - catalog[i].secfilt[Nsecfilt*j+Nsec].dM = stats.error; - catalog[i].secfilt[Nsecfilt*j+Nsec].Mstdev = 1000.0*stats.sigma; // Mstdev is in millimags (not enough space for more precision) - catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT; - catalog[i].secfilt[Nsecfilt*j+Nsec].Ncode = Ncode; - catalog[i].secfilt[Nsecfilt*j+Nsec].Nused = stats.Nmeas; - - catalog[i].secfilt[Nsecfilt*j+Nsec].M_80 = 1000 * stats.Upper80; - catalog[i].secfilt[Nsecfilt*j+Nsec].M_20 = 1000 * stats.Lower20; - catalog[i].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist; - - if ((Next > 0) && (Next > 0.5*N)) { - catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_OBJ_EXT; - } - } - } - } - - free (list); - free (dlist); - - free (aplist); - free (daplist); - return (TRUE); -} +# if (0) /* grab Nsec for named photcode */ @@ -776,4 +754,5 @@ return (TRUE); } +# endif /* set measure.Mcal for all measures except ID_MEAS_NOCAL and ID_IMAGE_PHOTOM_NOCAL */ @@ -800,5 +779,5 @@ if (isnan(Mgrid)) continue; - // note that measurements for which the image is not selected will not up modified + // note that measurements for which the image is not selected will not be modified // (that is a good thing! -- we keep a prior calibration) Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/initialize.c =================================================================== --- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/initialize.c (revision 33822) +++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/initialize.c (revision 33823) @@ -12,5 +12,21 @@ PhotcodeList = strcreate (argv[1]); photcodes = ParsePhotcodeList (PhotcodeList, &Nphotcodes, TRUE); // require SEC photcodes - } + } else { + 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); + } + } + if (USE_GRID && (Nphotcodes > 1)) { fprintf (stderr, "grid correction analysis currently can only operate on a single photcode\n"); Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/reload_catalogs.c =================================================================== --- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/reload_catalogs.c (revision 33822) +++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/reload_catalogs.c (revision 33823) @@ -86,5 +86,6 @@ TIMESTAMP(time4); - setMrelFinal (&catalog, flatcorr); + initMrel (&catalog, 1); + setMrelFinal (&catalog, flatcorr, FALSE); TIMESTAMP(time5); Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/relphot_objects.c =================================================================== --- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/relphot_objects.c (revision 33822) +++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/relphot_objects.c (revision 33823) @@ -90,8 +90,12 @@ } - setMave (&catalog, 1); + populate_tiny_values(&catalog, DVO_TV_MEASURE | DVO_TV_AVERAGE); + + initMrel (&catalog, 1); + setMrelFinal (&catalog, NULL, TRUE); if (!UPDATE) { dvo_catalog_unlock (&catalog); + free_tiny_values(&catalog); dvo_catalog_free (&catalog); continue; @@ -105,4 +109,5 @@ dvo_catalog_save (&catalog, VERBOSE); dvo_catalog_unlock (&catalog); + free_tiny_values(&catalog); dvo_catalog_free (&catalog); } Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/setMrelFinal.c =================================================================== --- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/setMrelFinal.c (revision 33822) +++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/setMrelFinal.c (revision 33823) @@ -5,5 +5,5 @@ // output dbFlags values -void setMrelFinal (Catalog *catalog, FlatCorrectionTable *flatcorr) { +void setMrelFinal (Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage) { off_t i; @@ -94,10 +94,11 @@ } - clean_measures (catalog, 1, TRUE, flatcorr); /* mark outliers ID_MEAS_POOR_PHOTOM */ + // 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 */ } - setMcalOutput (catalog, 1, flatcorr); + if (!simpleAverage) setMcalOutput (catalog, 1, flatcorr); /* clear ID_STAR_POOR, ID_STAR_FEW values before writing ??? */ @@ -258,5 +259,5 @@ } } - if (VERBOSE2) fprintf (stderr, "pass %d, Ntot: "OFF_T_FMT", Ntry: "OFF_T_FMT", Nskip: "OFF_T_FMT", Nkeep: "OFF_T_FMT"\n", + if (VERBOSE) fprintf (stderr, "pass %d, Ntot: "OFF_T_FMT", Ntry: "OFF_T_FMT", Nskip: "OFF_T_FMT", Nkeep: "OFF_T_FMT"\n", pass, Ntot, Ntry, Nskip, Nkeep); }