Index: /branches/eam_branches/ipp-20150419/Ohana/src/relphot/include/relphot.h =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/relphot/include/relphot.h (revision 38296) +++ /branches/eam_branches/ipp-20150419/Ohana/src/relphot/include/relphot.h (revision 38297) @@ -136,4 +136,5 @@ float *psfQfPerfMax; // max psfQfperf value for this secfilt + int *Nmeas; // count of PS1 exposure (chip) measurements 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 @@ -141,7 +142,4 @@ 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; @@ -396,4 +394,5 @@ void liststats_setmode PROTO((StatType *stats, char *strmode)); int liststats PROTO((double *value, double *dvalue, double *wvalue, int N, StatType *stats)); +int liststats_init PROTO((StatType *stats)); Catalog *load_catalogs PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath, char *syncfile)); Catalog *load_catalogs_parallel PROTO((SkyList *sky, int *Ncatalog, char *syncfile)); @@ -433,5 +432,5 @@ int setMrel PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); void setMrelFinal PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage)); -int setMrelOutput PROTO((Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr)); +int setMrelOutput PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); int setMave PROTO((Catalog *catalog, int Ncatalog)); void set_ZP PROTO((double ZERO)); @@ -495,8 +494,8 @@ 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 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); +int setMrel_catalog_alt (Catalog *catalog, int Nc, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt); +int setMrelAverageExposure (Catalog *catalog, int cat, off_t ave, int Nsecfilt, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results); +int setMrelAverageStack (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results); +int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results); int setGlobalObjStats (Average *average, Measure *measure); @@ -562,4 +561,7 @@ int whichGPC1filter (int photcode); int is2MASS (int photcode); +int isTYCHO (int photcode); + +int magStatsByRanking (StatDataSet *dataset, StatType *stats); SynthZeroPoints *SynthZeroPointsLoad (char *filename); @@ -570,2 +572,7 @@ int relphot_synthphot_catalog (Catalog *catalog, SynthZeroPoints *zpts); int relphot_synthphot_average (Average *average, SecFilt *secfilt, Measure *measure, SynthZeroPoints *zpts); + +void setMeasureRank (Catalog *catalog); +int getImageFlags (off_t meas, int cat); +int getMosaicFlags (off_t meas, int cat); + Index: /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/StarOps.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/StarOps.c (revision 38296) +++ /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/StarOps.c (revision 38297) @@ -111,4 +111,5 @@ results->psfQfPerfMax[i] = NAN; + results->Nmeas[i] = 0; results->Next[i] = 0; results->NexpPS1[i] = 0; @@ -121,6 +122,4 @@ results->kronData[i].Nlist = 0; } - results->NextPS1 = 0; - results->NpsfPS1 = 0; } @@ -149,4 +148,5 @@ ALLOCATE (results->psfQfPerfMax, float, results->Nsecfilt); + ALLOCATE (results->Nmeas, int, results->Nsecfilt); ALLOCATE (results->Next, int, results->Nsecfilt); ALLOCATE (results->NexpPS1, int, results->Nsecfilt); @@ -156,5 +156,5 @@ ALLOCATE (results->minUbercalDist, float, results->Nsecfilt); - SetMrelInfoResetObject (results, results->Nsecfilt); + SetMrelInfoResetObject (results); for (i = 0; i < NsynthNames; i++) { @@ -259,5 +259,5 @@ break; case 1: - setMrel_catalog_alt (catalog, i, -1, flatcorr, &results, Nsecfilt); + setMrel_catalog_alt (catalog, i, FALSE, flatcorr, &results, Nsecfilt); break; default: @@ -274,5 +274,5 @@ } -int setMrelOutput (Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr) { +int setMrelOutput (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { int i; @@ -293,5 +293,5 @@ // break; case 1: - setMrel_catalog_alt (catalog, i, pass, flatcorr, &results, Nsecfilt); + setMrel_catalog_alt (catalog, i, TRUE, flatcorr, &results, Nsecfilt); break; default: @@ -400,5 +400,5 @@ break; case 1: - setMrel_catalog_alt (catalog, i, -1, flatcorr, &results, Nsecfilt); + setMrel_catalog_alt (catalog, i, FALSE, flatcorr, &results, Nsecfilt); break; default: Index: /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/extra.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/extra.c (revision 38296) +++ /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/extra.c (revision 38297) @@ -77,2 +77,8 @@ } +int isTYCHO (int photcode) { + + if ((photcode == 2020) || (photcode == 2021)) return TRUE; + return FALSE; +} + Index: /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelCatalog.c (revision 38296) +++ /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelCatalog.c (revision 38297) @@ -1,4 +1,5 @@ # include "relphot.h" # define UBERCAL_WEIGHT 100.0 +void sort_StatDataSet (StatDataSet *dataset); # define SKIP_THIS_MEAS(REASON) { \ @@ -16,8 +17,10 @@ # define CHECK_VALID_MAG(MAG) (!isnan(MAG) && (MAG > 0.0) && (MAG < 30.0)) +# if (0) static float MagToFlux (float Mag) { float Flux = pow(10.0, -0.4*(Mag)); return (Flux); } +# endif int print_measure_set_alt (Average *average, SecFilt *secfilt, Measure *measure) { @@ -38,5 +41,5 @@ } -int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) { +int setMrel_catalog_alt (Catalog *catalog, int Nc, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) { off_t j; @@ -48,24 +51,16 @@ SetMrelInfoReset (results); // reset the countesrs - int isSetMrelFinal = (pass >= 0); - for (j = 0; j < catalog[Nc].Naverage; j++) { - // we are guaranteed to have averageT, but not average - off_t m = catalog[Nc].averageT[j].measureOffset; - Average *average = catalog[Nc].average ? &catalog[Nc].average[j] : NULL; - Measure *measure = catalog[Nc].measure ? &catalog[Nc].measure[m] : NULL; - SecFilt *secfilt = catalog[Nc].secfilt ? &catalog[Nc].secfilt[j*Nsecfilt] : NULL; - setMrelAverageExposure (m, Nc, isSetMrelFinal, flatcorr, results, average, &catalog[Nc].averageT[j], secfilt, measure, &catalog[Nc].measureT[m], &catalog[Nc].found_t[Nsecfilt*j]); + + setMrelAverageExposure (&catalog[Nc], Nc, j, Nsecfilt, isSetMrelFinal, flatcorr, results); // only apply Stack operation on setMrelFinal in first pass - if (isSetMrelFinal && (pass == 0) && !IS_DIFF_DB) { - setMrelAverageStack (m, Nc, flatcorr, results, average, secfilt, &catalog[Nc].measure[m], catalog[Nc].Nmeasure); - setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]); + if (isSetMrelFinal && !IS_DIFF_DB) { + setMrelAverageStack (&catalog[Nc], Nc, j, Nsecfilt, flatcorr, results); } // only measure force-warp mean values if issetMrelFinal (make it optional?) if (isSetMrelFinal) { - setMrelAverageForcedWarp (m, Nc, pass, results, average, secfilt, measure, &catalog[Nc].foundWarp_t[Nsecfilt*j], catalog[Nc].Nmeasure); - setGlobalObjStats (average, measure); + setMrelAverageForcedWarp (&catalog[Nc], Nc, j, Nsecfilt, flatcorr, results); } } @@ -91,8 +86,23 @@ // 4) some reference photcode of some kind can be specified as fixed and have a high weight +static SynthZeroPoints *synthzpts = NULL; + // set mean of chip measurements (selected by photcode range for now): -int setMrelAverageExposure (off_t measureOffset, int cat, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT) { - +int setMrelAverageExposure (Catalog *catalog, int cat, off_t ave, int Nsecfilt, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results) { + + int Nsec; off_t k; + + // we are guaranteed to have averageT, but not average + AverageTiny *averageT = &catalog[cat].averageT[ave]; + + off_t measureOffset = averageT->measureOffset; + MeasureTiny *measureT = &catalog[cat].measureT[measureOffset]; + + // we are NOT guaranteed to have average, measure, secfilt + Average *average = catalog[cat].average ? &catalog[cat].average[ave] : NULL; + Measure *measure = catalog[cat].measure ? &catalog[cat].measure[measureOffset] : NULL; + SecFilt *secfilt = catalog[cat].secfilt ? &catalog[cat].secfilt[ave*Nsecfilt] : NULL; + char *measureRank = catalog[cat].measureRank ? &catalog[cat].measureRank[measureOffset] : NULL; // we are measuring means for 3 types of mags: psf, ap, kron. @@ -110,4 +120,11 @@ SetMrelInfoResetObject (results); // reset the per-object arrays + + int NextPS1 = 0; + int NpsfPS1 = 0; + int GoodPS1 = FALSE; + int Good2MASS = FALSE; + int Galaxy2MASS = FALSE; + int haveTYCHO = FALSE; // assign measurements to the photcode lists @@ -124,5 +141,5 @@ if (isGPC1warp(measureT[k].photcode)) continue; - if (isTycho(measureT[k].photcode)) { haveTYCHO = TRUE; } + if (isTYCHO(measureT[k].photcode)) { haveTYCHO = TRUE; } if (is2MASS(measureT[k].photcode)) { @@ -138,4 +155,6 @@ int Nsec = GetPhotcodeNsec (code->equiv); if (Nsec < 0) continue; // skip measurements which do not have an equiv average photcode + + results->Nmeas[Nsec] ++; // force the use of SYN even if we have PS1 mags? @@ -148,7 +167,9 @@ // 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; - + if (isSetMrelFinal) { + if (measure[k].psfQF > results->psfQfMax[Nsec]) results->psfQfMax[Nsec] = measure[k].psfQF; + if (measure[k].psfQFperf > results->psfQfPerfMax[Nsec]) results->psfQfPerfMax[Nsec] = measure[k].psfQFperf; + } + // are we a PS1 exposure photcode? (hard-wired photcodes) if (isGPC1chip(measureT[k].photcode)) { results->NexpPS1[Nsec] ++; results->havePS1[Nsec] = TRUE; } @@ -156,5 +177,6 @@ // ** 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) { + off_t meas = measureOffset + k; + 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 @@ -173,5 +195,5 @@ // XXX ubercal dist should be by Nsec int myUbercalDist = getUbercalDist(meas, cat); - result->minUbercalDist[Nsec] = MIN(result->minUbercalDist[Nsec], myUbercalDist); + results->minUbercalDist[Nsec] = MIN(results->minUbercalDist[Nsec], myUbercalDist); int isUbercal = (measureT[k].dbFlags & ID_MEAS_PHOTOM_UBERCAL); @@ -184,14 +206,13 @@ if (isSetMrelFinal) { Map = PhotCat (&measure[k], MAG_CLASS_APER); - dMap = MAX (hypot(measureT[k].dMap, code->photomErrSys), MIN_ERROR); + dMap = MAX (hypot(measure[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 ++; + int Nap = results->aperData[Nsec].Nlist; + results->aperData[Nsec].flxlist[Nap] = Map - Mcal - Mmos - Mgrid; + results->aperData[Nsec].errlist[Nap] = dMap; + results->aperData[Nsec].wgtlist[Nap] = useUbercalWeight ? UBERCAL_WEIGHT : 1.0; + results->aperData[Nsec].ranking[Nap] = measureRank[k]; + results->aperData[Nsec].measSeq[Nap] = k; + results->aperData[Nsec].Nlist ++; } } @@ -200,9 +221,9 @@ if (isSetMrelFinal) { Mkron = PhotCat (&measure[k], MAG_CLASS_KRON); - dMkron = MAX (hypot(measureT[k].dMkron, code->photomErrSys), MIN_ERROR); + dMkron = MAX (hypot(measure[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].errlist[Nkron] = dMkron; results->kronData[Nsec].wgtlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0; results->kronData[Nsec].ranking[Nkron] = measureRank[k]; @@ -237,7 +258,7 @@ if (dMagAp > SigmaAp) { results->Next[Nsec] ++; - results->NextPS1 ++; + NextPS1 ++; } else { - results->NpsfPS1 ++; + NpsfPS1 ++; } } @@ -246,11 +267,15 @@ // now calculate the mean stats for the Nsec bands. for (Nsec = 0; Nsec < Nsecfilt; Nsec++) { - dvo_init_secfilt (&secfilt[Nsec]); + dvo_secfilt_init (&secfilt[Nsec]); // XXX does this reset astrometry bits? no, right? // XXX hardwired - if ((thisCode < 6) || (thisCode == 9)) { + if ((Nsec < 4) || (Nsec == 8)) { secfilt[Nsec].Ncode = results->NexpPS1[Nsec]; } else { - secfilt[Nsec].Ncode = results->Ncode[Nsec]; // 2MASS data if it exists + secfilt[Nsec].Ncode = results->Nmeas[Nsec]; // 2MASS data if it exists + } + + if (haveTYCHO) { + // secfilt[Nsec].flags |= ID_SECF_HAVE_TYCHO; } @@ -263,25 +288,25 @@ float ZP = 0.0; - if (!haveTycho) { + // if we have loaded the synth zpt correction, and do NOT have a tycho measurement, apply the correction + if (!haveTYCHO && synthzpts) { // 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); + RD_to_XY (&X, &Y, R, measure->D, &synthzpts->coords); if (X < 0) return FALSE; if (Y < 0) return FALSE; - if (X >= zpts->Nx) return FALSE; - if (Y >= zpts->Ny) return FALSE; + if (X >= synthzpts->Nx) return FALSE; + if (Y >= synthzpts->Ny) return FALSE; int Xpix = X; int Ypix = Y; - int Npix = Xpix + Ypix*zpts->Nx; + int Npix = Xpix + Ypix*synthzpts->Nx; - float *value = (float *) zpts->matrix[Nsec].buffer; + float *value = (float *) synthzpts->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].Mchisq = 0.0; secfilt[Nsec].flags |= ID_SECF_USE_SYNTH; continue; @@ -289,5 +314,6 @@ // if too few valid measurements meet the minimum criteria, go to the next entry - int Nranking = magStatsByRanking (&results->psfData[Nsec], &results->psfstats); + StatType *psfstats = &results->psfstats; + int Nranking = magStatsByRanking (&results->psfData[Nsec], psfstats); if (Nranking < Nminmeas) { secfilt[Nsec].flags |= ID_STAR_FEW; @@ -300,31 +326,39 @@ // when running -averages, we have no information about the images, so we cannot set this - if (minUbercalDist > -1) { - secfilt[Nsec].ubercalDist = minUbercalDist; + if (results->minUbercalDist[Nsec] > -1) { + secfilt[Nsec].ubercalDist = results->minUbercalDist[Nsec]; } 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; - + if (Nranking) { + secfilt[Nsec].Mstdev = psfstats->sigma; + secfilt[Nsec].Nused = psfstats->Nmeas; + secfilt[Nsec].Mmax = psfstats->max; + secfilt[Nsec].Mmin = psfstats->min; + } + + secfilt[Nsec].psfQfMax = results->psfQfMax[Nsec]; + secfilt[Nsec].psfQfPerfMax = results->psfQfPerfMax[Nsec]; + + StatType *apstats = &results->apstats; + Nranking = magStatsByRanking (&results->aperData[Nsec], apstats); + if (Nranking) { + secfilt[Nsec].Map = apstats->mean; + secfilt[Nsec].dMap = apstats->error; + secfilt[Nsec].sMap = apstats->sigma; + secfilt[Nsec].NusedAp = Nranking; + } + + StatType *kronstats = &results->kronstats; 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; + if (Nranking) { + 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])) { + if ((results->Next[Nsec] > 0) && (results->Next[Nsec] > 0.5*psfstats->Nmeas)) { secfilt[Nsec].flags |= ID_SECF_OBJ_EXT; } @@ -362,5 +396,5 @@ if (isSetMrelFinal) { - if (results->NextPS1 && (results->NextPS1 > results->NpsfPS1)) { + if (NextPS1 && (NextPS1 > NpsfPS1)) { average[0].flags |= ID_OBJ_EXT; } @@ -378,4 +412,5 @@ } +# if (0) // 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) { @@ -750,4 +785,5 @@ return (TRUE); } +# endif // only apply Stack operation on setMrelFinal in first pass @@ -756,9 +792,17 @@ // 2) select the BEST detections per filter (regardless of PRIMARY) // 3) apply the zero point and AB->Jy transformations -int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t Nmeasure) { +int setMrelAverageStack (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results) { + + // we are guaranteed to have average, measure, secfilt + Average *average = &catalog[cat].average[ave]; + off_t measureOffset = average->measureOffset; + int Nmeasure = average->Nmeasure; + + Measure *measure = &catalog[cat].measure[measureOffset]; + SecFilt *secfilt = &catalog[cat].secfilt[ave*Nsecfilt]; off_t k; - float Mcal= 0, Mmos = 0, Mgrid = 0, Finst = 0; + float Mcal = 0, Mmos = 0, Mgrid = 0, Finst = 0; // set the primary projection cell and skycell for this coordinate @@ -801,5 +845,5 @@ off_t meas = measureOffset; - for (k = 0; k < average[0].Nmeasure; k++, meas++) { + for (k = 0; k < Nmeasure; k++, meas++) { // only examine gpc1 stack data @@ -962,4 +1006,7 @@ average[0].flags |= ID_OBJ_BAD_STACK; } + + setGlobalObjStats (average, measure); + return (TRUE); } @@ -974,8 +1021,17 @@ // somewhat simplified relative to chip-photometry: // * no grid, no mosaic, no 2MASS, no SYNTH, no Ubercal, no flatcorr -int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found, off_t Nmeasure) { +int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results) { + +# if (0) + + off_t measureOffset = average->measureOffset; + + // we are guaranteed to have average, measure, secfilt + Average *average = &catalog[cat].average[ave]; + Measure *measure = &catalog[cat].measure[measureOffset]; + SecFilt *secfilt = &catalog[cat].secfilt[ave*Nsecfilt]; off_t k; - float Mcal= 0; + float Mcal = 0; // we are measuring means for 3 types of FLUXes: psf, ap, kron. I am using the psf mag @@ -1180,4 +1236,9 @@ secfilt[Nsec].NwarpGood = NwarpGood; } + + setGlobalObjStats (average, measure); + +# endif + return (TRUE); } @@ -1220,5 +1281,5 @@ // 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); + sort_StatDataSet (dataset); // find the values with the same minimum rank: @@ -1229,10 +1290,25 @@ 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); } + +/* sort a coordinate pair (X,Y) and the associated index (S) */ +void sort_StatDataSet (StatDataSet *dataset) { + +# define SWAPFUNC(A,B){ double dtmp; int itmp; \ + dtmp = dataset->flxlist[A]; dataset->flxlist[A] = dataset->flxlist[B]; dataset->flxlist[B] = dtmp; \ + dtmp = dataset->errlist[A]; dataset->errlist[A] = dataset->errlist[B]; dataset->errlist[B] = dtmp; \ + dtmp = dataset->wgtlist[A]; dataset->wgtlist[A] = dataset->wgtlist[B]; dataset->wgtlist[B] = dtmp; \ + itmp = dataset->ranking[A]; dataset->ranking[A] = dataset->ranking[B]; dataset->ranking[B] = itmp; \ + itmp = dataset->measSeq[A]; dataset->measSeq[A] = dataset->measSeq[B]; dataset->measSeq[B] = itmp; \ + } +# define COMPARE(A,B)(dataset->ranking[A] < dataset->ranking[B]) + + OHANA_SORT (dataset->Nlist, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} + Index: /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelFinal.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelFinal.c (revision 38296) +++ /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelFinal.c (revision 38297) @@ -102,6 +102,6 @@ ALLOCATE (catalog[0].measureRank, char, catalog[0].Nmeasure); - setMeasurementRank (catalog); - setMrelOutput (catalog); + setMeasureRank (catalog); + setMrelOutput (catalog, 1, flatcorr); if (!simpleAverage) setMcalOutput (catalog, 1, flatcorr); @@ -311,5 +311,5 @@ // measurements ranked by inst mag limit (not REF mags) if (ImagSelect) { - mag = PhotInst (&measure[i], MAG_CLASS_PSF); + float mag = PhotInst (&measure[i], MAG_CLASS_PSF); if (mag < ImagMin) { measureRank[i] = 8; continue; } if (mag > ImagMax) { measureRank[i] = 8; continue; }