Index: /branches/eam_branches/ipp-20140206/Ohana/src/relphot/include/relphot.h =================================================================== --- /branches/eam_branches/ipp-20140206/Ohana/src/relphot/include/relphot.h (revision 36534) +++ /branches/eam_branches/ipp-20140206/Ohana/src/relphot/include/relphot.h (revision 36535) @@ -95,10 +95,15 @@ // 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 *Mpsflist; // list of measure.mag values for a given star + double *dpsflist; // mag errors for a star + double *wpsflist; // weights to use for mean mags + + double *Maplist; // ap mags for a star + double *daplist; // ap mags for a star + double *waplist; // ap mags for a star + + double *Mkronlist; // kron mags for a star double *dkronlist; // kron mag errors + double *wkronlist; // kron mag errors double *psfqf_list; // psfqf for all filters @@ -106,5 +111,5 @@ double *stargal_list; // stargal for all filters - StatType stats; + StatType psfstats; StatType apstats; StatType kronstats; Index: /branches/eam_branches/ipp-20140206/Ohana/src/relphot/src/StarOps.c =================================================================== --- /branches/eam_branches/ipp-20140206/Ohana/src/relphot/src/StarOps.c (revision 36534) +++ /branches/eam_branches/ipp-20140206/Ohana/src/relphot/src/StarOps.c (revision 36535) @@ -63,19 +63,19 @@ results->Ngrid = 0; if (allocLists) { - ALLOCATE (results->list, double, Nmax); - ALLOCATE (results->dlist, double, Nmax); - ALLOCATE (results->wlist, double, Nmax); + ALLOCATE (results->Mpsflist, double, Nmax); + ALLOCATE (results->dpsflist, double, Nmax); + ALLOCATE (results->wpsflist, double, Nmax); } } void SetMrelInfoFree (SetMrelInfo *results) { - free (results->list); - free (results->dlist); - free (results->wlist); + free (results->Mpsflist); + free (results->dpsflist); + free (results->wpsflist); } void SetMrelInfoAccum (SetMrelInfo *summary, SetMrelInfo *results) { summary->Nfew += results->Nfew ; - summary->Ncode += results->Ncode ; + summary->Ncode += results->Ncode ; summary->Nsys += results->Nsys ; summary->Nbad += results->Nbad ; @@ -158,7 +158,12 @@ SetMrelInfoInit (&summary, FALSE); SetMrelInfoInit (&results, TRUE); // allocates results->list,dlist,wlist - ALLOCATE (results.aplist, double, Nmax); - ALLOCATE (results.kronlist, double, Nmax); + + ALLOCATE (results.Maplist, double, Nmax); + ALLOCATE (results.daplist, double, Nmax); + ALLOCATE (results.waplist, double, Nmax); + + ALLOCATE (results.Mkronlist, double, Nmax); ALLOCATE (results.dkronlist, double, Nmax); + ALLOCATE (results.wkronlist, double, Nmax); ALLOCATE (results.psfqf_list, double, Nmax); @@ -183,7 +188,11 @@ SetMrelInfoFree (&results); - free (results.aplist); - free (results.kronlist); + free (results.Maplist); + free (results.daplist); + free (results.waplist); + + free (results.Mkronlist); free (results.dkronlist); + free (results.wkronlist); free (results.psfqf_list); @@ -331,4 +340,6 @@ exit (3); +# if (0) + off_t j, k, m, ID; int N; @@ -336,5 +347,5 @@ StatType stats, apstats, kronstats; - liststats_setmode (&stats, STATMODE); + liststats_setmode (&psfstats, STATMODE); liststats_setmode (&apstats, STATMODE); liststats_setmode (&kronstats, STATMODE); @@ -837,4 +848,5 @@ if (primaryCell) free (primaryCell); return (TRUE); +# endif } Index: /branches/eam_branches/ipp-20140206/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- /branches/eam_branches/ipp-20140206/Ohana/src/relphot/src/setMrelCatalog.c (revision 36534) +++ /branches/eam_branches/ipp-20140206/Ohana/src/relphot/src/setMrelCatalog.c (revision 36535) @@ -34,5 +34,5 @@ off_t j; - liststats_setmode (&results->stats, STATMODE); + liststats_setmode (&results->psfstats, STATMODE); liststats_setmode (&results->apstats, STATMODE); liststats_setmode (&results->kronstats, STATMODE); @@ -67,16 +67,23 @@ 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; + // 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; - - StatType *stats = &results->stats; - StatType *apstats = &results->apstats; + double *wkronlist = results->wkronlist; + + StatType *psfstats = &results->psfstats; + StatType *apstats = &results->apstats; StatType *kronstats = &results->kronstats; @@ -125,5 +132,7 @@ off_t meas = measureOffset; - N = 0; + int Nap = 0; + int Npsf = 0; + int Nkron = 0; for (k = 0; k < averageT[0].Nmeasure; k++, meas++) { @@ -164,4 +173,9 @@ 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) { @@ -231,5 +245,5 @@ if (measureT[k].M < MaxMagForceSynth) { forceSynth = TRUE; - forceSynthEntry = N; + forceSynthEntry = Npsf; } else { if (pass < 4) { @@ -240,12 +254,21 @@ } - // Map (aplist) and Mkron (kronlist,dkronlist) are used to calculate mean mags per filter + // Map (Maplist) and Mkron (Mkronlist,dkronlist) are used to calculate mean mags per filter float Map = PhotAper (&measure[k]); - aplist[N] = Map - Mcal - Mmos - Mgrid; + if (!isnan(Map)) { + Maplist[Nap] = Map - Mcal - Mmos - Mgrid; + daplist[Nap] = dMpsf; + waplist[Nap] = isUbercal ? 10.0 : 1.0; + Nap ++; + } float Mkron = PhotKron (&measure[k]); - kronlist[N] = Mkron - Mcal - Mmos - Mgrid; - dkronlist[N] = measure[k].dMkron; - } + if (!isnan(Mkron)) { + Mkronlist[Nkron] = Mkron - Mcal - Mmos - Mgrid; + dkronlist[Nkron] = measure[k].dMkron; + wkronlist[Nkron] = isUbercal ? 10.0 : 1.0; + Nkron ++; + } + } // if (isSetMrelFinal) // dlist gives the error per measurement, wlist gives the weight @@ -255,7 +278,7 @@ // 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; + Mpsflist[Npsf] = Msys - Mcal - Mmos - Mgrid; + dpsflist[Npsf] = dMpsf; + wpsflist[Npsf] = isUbercal ? 10.0 : 1.0; // NOTE: @@ -264,19 +287,14 @@ // Mmos and Mgrid are offsets for mosaic and grid - // 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 + // XXX probably need to use the photocde table to assign reference mag weights. if (refPhotcode) { if (code->code == refPhotcode->code) { - wlist[N] = 100.0; - } - } - N++; + wpsflist[Npsf] = 100.0; + } + } + Npsf ++; } @@ -293,6 +311,6 @@ // 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); + 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) { @@ -308,15 +326,15 @@ 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; + 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].Xm = (psfstats->Nmeas > 1) ? 100.0*log10(psfstats->chisq + 1e-4) : NAN_S_SHORT; // when running -averages, we have no information about the images, so we cannot set this @@ -328,21 +346,21 @@ found[Nsec] = TRUE; - secfilt[Nsec].Mstdev = 1000.0*stats->sigma; // Mstdev is in millimags (not enough space for more precision) + secfilt[Nsec].Mstdev = 1000.0*psfstats->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; + secfilt[Nsec].Nused = psfstats->Nmeas; + + secfilt[Nsec].M_80 = 1000 * psfstats->Upper80; + secfilt[Nsec].M_20 = 1000 * psfstats->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; + liststats (Maplist, daplist, waplist, Nap, apstats); + secfilt[Nsec].Map = Nap > 0 ? apstats->mean : NAN; + + liststats (Mkronlist, dkronlist, wkronlist, Nkron, kronstats); + secfilt[Nsec].Mkron = Nkron > 0 ? kronstats->mean : NAN; + secfilt[Nsec].dMkron = Nkron > 0 ? kronstats->error : NAN; // NOTE: for 2MASS measurements, Next should be 1, as should N - if ((Next > 0) && (Next > 0.5*N)) { + if ((Next > 0) && (Next > 0.5*Npsf)) { secfilt[Nsec].flags |= ID_SECF_OBJ_EXT; }