Index: /branches/eam_branches/ohana.20170822/src/relphot/include/relphot.h =================================================================== --- /branches/eam_branches/ohana.20170822/src/relphot/include/relphot.h (revision 40212) +++ /branches/eam_branches/ohana.20170822/src/relphot/include/relphot.h (revision 40213) @@ -425,6 +425,5 @@ 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 getMcal PROTO((off_t meas, int cat, dvoMagClassType class)); float getMflat PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog)); float getMgrid PROTO((off_t meas, int cat)); Index: /branches/eam_branches/ohana.20170822/src/relphot/src/GridOps.c =================================================================== --- /branches/eam_branches/ohana.20170822/src/relphot/src/GridOps.c (revision 40212) +++ /branches/eam_branches/ohana.20170822/src/relphot/src/GridOps.c (revision 40213) @@ -401,5 +401,5 @@ // skip images marked as BAD - Mcal = getMcal (m, c, flatcorr, catalog); + Mcal = getMcal (m, c, MAG_CLASS_PSF); if (isnan(Mcal)) { Ncal ++; @@ -558,5 +558,5 @@ continue; } - Mcal = getMcal (m, c, flatcorr, catalog); + Mcal = getMcal (m, c, MAG_CLASS_PSF); if (isnan(Mcal)) { Ncal ++; @@ -645,5 +645,5 @@ continue; } - Mcal = getMcal (m, c, flatcorr, catalog); + Mcal = getMcal (m, c, MAG_CLASS_PSF); if (isnan(Mcal)) continue; Mmos = getMmos (m, c); Index: /branches/eam_branches/ohana.20170822/src/relphot/src/ImageOps.c =================================================================== --- /branches/eam_branches/ohana.20170822/src/relphot/src/ImageOps.c (revision 40212) +++ /branches/eam_branches/ohana.20170822/src/relphot/src/ImageOps.c (revision 40213) @@ -432,55 +432,25 @@ } -// returns image.Mcal +// returns image.McalPSF or image.McalAPER // NOTE: static flat-field component is included in measure.Mflat -float getMcal (off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog) { - OHANA_UNUSED_PARAM(flatcorr); - OHANA_UNUSED_PARAM(catalog); - - off_t i; - float value; +float getMcal (off_t meas, int cat, dvoMagClassType class) { + + off_t i; i = MeasureToImage[cat][meas]; - if (i == -1) return (NAN); - - if (image[i].flags & IMAGE_BAD) return (NAN); - value = image[i].McalPSF; - - // to do this, I need to pass in the catalog and flatcorr pointers - // int flat_id = image[i].photom_map_id; - // float offset = 0.0; - // 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) -// NOTE: static flat-field component is included in measure.Mflat -float getMcal_alt (off_t meas, int cat, FlatCorrectionTable *flatcorr, float Xccd, float Yccd) { - OHANA_UNUSED_PARAM(flatcorr); - OHANA_UNUSED_PARAM(Xccd); - OHANA_UNUSED_PARAM(Yccd); - - off_t i; - float value; - - i = MeasureToImage[cat][meas]; - if (i == -1) return (NAN); - - if (image[i].flags & IMAGE_BAD) return (NAN); - value = image[i].McalPSF; - - // to do this, I need to pass in the catalog and flatcorr pointers - // int flat_id = image[i].photom_map_id; - // float offset = 0.0; - // if (flat_id) { - // offset = FlatCorrectionOffset (flatcorr, flat_id, Xccd, Yccd); - // } - // value -= offset; - - return (value); + if (i == -1) return NAN; + + if (image[i].flags & IMAGE_BAD) return NAN; + + switch (class) { + case MAG_CLASS_PSF: + return image[i].McalPSF; + case MAG_CLASS_APER: + case MAG_CLASS_KRON: + return image[i].McalAPER; + default: + return NAN; + } + return NAN; // should not be able to reach here } @@ -615,8 +585,6 @@ void setMcal (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { - off_t i, j, m, c, n, Nmax; + off_t i, j, m, c, n; int mark, bad, Nfew, Nbad, Nmos, Nrel, Ngrid, Nsys; - float Msys, Mrel, Mmos, Mgrid, Mflat; - // double *list, *dlist, *Mlist, *dMlist; StatType stats; @@ -635,9 +603,11 @@ } - Nmax = 0; + off_t Nmax = 0; for (i = 0; i < Nimage; i++) { Nmax = MAX (Nmax, N_onImage[i]); } - StatDataSet *refStars = StatDataSetAlloc (1, Nmax); + + StatDataSet *kronStars = StatDataSetAlloc (1, Nmax); + StatDataSet *psfStars = StatDataSetAlloc (1, Nmax); StatDataSet *brightStars = StatDataSetAlloc (1, Nmax); @@ -690,10 +660,10 @@ continue; } - Mmos = getMmos (m, c); + float Mmos = getMmos (m, c); if (isnan(Mmos)) { Nmos ++; continue; } - Mgrid = getMgrid (m, c); + float Mgrid = getMgrid (m, c); if (isnan(Mgrid)) { Ngrid++; @@ -701,5 +671,5 @@ } - // Mrel is the average magnitude for this star. For PS1 stacks, we have too much + // Mrel* is the average magnitude for this star. For PS1 stacks, we have too much // PSF variability. We need to calibrate the PSF magnitudes separately from the // Aperture-like magnitues. (We have an option to use the kron magnitudes or the @@ -707,9 +677,10 @@ // magnitude type - Mrel = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); - if (isnan(Mrel)) { + float MrelPSF = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); + if (isnan(MrelPSF)) { Nrel ++; continue; } + float MrelKron = getMrel (catalog, m, c, MAG_CLASS_KRON, MAG_SRC_CHP); // image.Mcal is not supposed to include the flat-field correction, so we need to @@ -718,12 +689,13 @@ // the flat-correction. NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat) - Mflat = getMflat (m, c, flatcorr, catalog); + float Mflat = getMflat (m, c, flatcorr, catalog); n = catalog[c].measureT[m].averef; - Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); - if (isnan(Msys)) { + float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); + if (isnan(MsysPSF)) { Nsys++; continue; } + float MsysKron = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_KRON); PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode); @@ -735,19 +707,24 @@ skip: - refStars->flxlist[Nref] = Msys - Mrel - Mmos - Mgrid + Mflat; - refStars->errlist[Nref] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); - refStars->wgtlist[Nref] = 1; - refStars->msklist[Nref] = 0; - if (fabs(refStars->flxlist[Nref]) > 0.03) { + psfStars->flxlist[Nref] = MsysPSF - MrelPSF - Mmos - Mgrid + Mflat; + psfStars->errlist[Nref] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); + psfStars->wgtlist[Nref] = 1; + psfStars->msklist[Nref] = 0; + if (fabs(psfStars->flxlist[Nref]) > 0.03) { // fprintf (stderr, "deviant\n"); } + kronStars->flxlist[Nref] = MsysKron - MrelKron - Mmos - Mgrid + Mflat; + kronStars->errlist[Nref] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); + kronStars->wgtlist[Nref] = 1; + kronStars->msklist[Nref] = 0; + if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) { - fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f %6.3f %6.3f %6.3f %6.3f : %6.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, Msys, Mrel, Mmos, Mgrid, Mflat, refStars->flxlist[Nref]); + fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f %6.3f %6.3f %6.3f %6.3f : %6.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, Mmos, Mgrid, Mflat, psfStars->flxlist[Nref]); } if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) { - brightStars->flxlist[Nbright] = refStars->flxlist[Nref]; - brightStars->errlist[Nbright] = refStars->errlist[Nref]; + brightStars->flxlist[Nbright] = psfStars->flxlist[Nref]; + brightStars->errlist[Nbright] = psfStars->errlist[Nref]; brightStars->wgtlist[Nbright] = 1; brightStars->msklist[Nbright] = 0; @@ -774,10 +751,9 @@ // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not) # if (BASIC_STATS) - liststats (refStars->flxlist, refStars->errlist, NULL, Nref, &stats); + liststats (psfStars->flxlist, psfStars->errlist, NULL, Nref, &stats); # else - liststats_irls (refStars, Nref, &stats); + liststats_irls (psfStars, Nref, &stats); # endif image[i].McalPSF = stats.mean; - image[i].McalAPER = stats.mean; image[i].dMcal = stats.error; image[i].nFitPhotom = Nref; @@ -785,7 +761,15 @@ Ncalibrated ++; + // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not) +# if (BASIC_STATS) + liststats (kronStars->flxlist, kronStars->errlist, NULL, Nref, &stats); +# else + liststats_irls (kronStars, Nref, &stats); +# endif + image[i].McalAPER = stats.mean; + if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) { for (j = 0; j < Nref; j++) { - fprintf (stderr, "%1d, %8d : %6.3f %6.3f %6.3f %d\n", (int) i, (int) image[i].imageID, refStars->flxlist[j], refStars->errlist[j], refStars->wgtlist[j], refStars->msklist[j]); + fprintf (stderr, "%1d, %8d : %6.3f %6.3f %6.3f %d\n", (int) i, (int) image[i].imageID, psfStars->flxlist[j], psfStars->errlist[j], psfStars->wgtlist[j], psfStars->msklist[j]); } } @@ -801,5 +785,5 @@ if (PLOTSTUFF) { fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].McalPSF, image[i].dMcal); - plot_setMcal (refStars->flxlist, Nref, &stats, CLOUD_TOLERANCE); + plot_setMcal (psfStars->flxlist, Nref, &stats, CLOUD_TOLERANCE); } @@ -818,5 +802,5 @@ StatDataSetFree (brightStars, 1); - StatDataSetFree (refStars, 1); + StatDataSetFree (psfStars, 1); fprintf (stderr, "%d images calibrated\n", Ncalibrated); @@ -1056,5 +1040,5 @@ c = ImageToCatalog[i][j]; - Mcal = getMcal (m, c, flatcorr, catalog); + Mcal = getMcal (m, c, MAG_CLASS_PSF); if (isnan(Mcal)) continue; Mmos = getMmos (m, c); Index: /branches/eam_branches/ohana.20170822/src/relphot/src/MosaicOps.c =================================================================== --- /branches/eam_branches/ohana.20170822/src/relphot/src/MosaicOps.c (revision 40212) +++ /branches/eam_branches/ohana.20170822/src/relphot/src/MosaicOps.c (revision 40213) @@ -925,8 +925,10 @@ off_t Nmax; int PoorImages; - double *list; - double *dlist; - double *Mlist; - double *dMlist; + double *psfMagList; + double *psfErrList; + double *brightMagList; + double *brightErrList; + double *kronMagList; + double *kronErrList; } SetMmosInfo; @@ -959,16 +961,20 @@ if (allocLists) { - ALLOCATE (info->list, double, Nmax); - ALLOCATE (info->dlist, double, Nmax); - ALLOCATE (info->Mlist, double, Nmax); - ALLOCATE (info->dMlist, double, Nmax); + ALLOCATE (info->psfMagList, double, Nmax); + ALLOCATE (info->psfErrList, double, Nmax); + ALLOCATE (info->kronMagList, double, Nmax); + ALLOCATE (info->kronErrList, double, Nmax); + ALLOCATE (info->brightMagList, double, Nmax); + ALLOCATE (info->brightErrList, double, Nmax); } } void SetMmosInfoFree (SetMmosInfo *info) { - free (info->list); - free (info->dlist); - free (info->Mlist); - free (info->dMlist); + free (info->psfMagList ); + free (info->psfErrList ); + free (info->kronMagList ); + free (info->kronErrList ); + free (info->brightMagList); + free (info->brightErrList); } @@ -1071,8 +1077,10 @@ liststats_setmode (&stats, "INNER_WTMEAN"); - double *list = info->list; - double *dlist = info->dlist; - double *Mlist = info->Mlist; - double *dMlist = info->dMlist; + double *psfMagList = info->psfMagList; + double *psfErrList = info->psfErrList; + double *kronMagList = info->kronMagList; + double *kronErrList = info->kronErrList; + double *brightMagList = info->brightMagList; + double *brightErrList = info->brightErrList; assert (Nmos >= 0); @@ -1139,5 +1147,4 @@ int N = 0; for (j = 0; j < N_onMosaic[Nmos]; j++) { - float Msys, Mrel, Mcal, Mgrid, Mflat; off_t m = MosaicToMeasure[Nmos][j]; @@ -1146,17 +1153,17 @@ // NOTE : we are only using Mcal == McalPSF in this function; we set McalAPER to McalPSF if (fout) { - Mcal = getMcal (m, c, flatcorr, catalog); - Mgrid = getMgrid (m, c); - Mrel = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); - Mflat = getMflat (m, c, flatcorr, catalog); + float Mcal = getMcal (m, c, MAG_CLASS_PSF); + float Mgrid = getMgrid (m, c); + float MrelPSF = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); + float Mflat = getMflat (m, c, flatcorr, catalog); off_t n = catalog[c].measureT[m].averef; - Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); - - float delta = Msys - Mrel - Mcal - Mgrid + Mflat; + float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); + + float delta = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat; int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD); - fprintf (fout, "%f %f : %f %f %f %f %f : %f %d\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, Msys, Mrel, Mcal, Mgrid, Mflat, delta, isBad); + fprintf (fout, "%f %f : %f %f %f %f %f : %f %d\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, Mcal, Mgrid, Mflat, delta, isBad); } @@ -1165,19 +1172,27 @@ continue; } - Mcal = getMcal (m, c, flatcorr, catalog); + float Mcal = getMcal (m, c, MAG_CLASS_PSF); if (isnan(Mcal)) { info->Ncal++; continue; } - Mgrid = getMgrid (m, c); + float Mgrid = getMgrid (m, c); if (isnan(Mgrid)) { info->Ngrid ++; continue; } - Mrel = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); - if (isnan(Mrel)) { + + // Mrel* is the average magnitude for this star. For PS1 stacks, we have too much + // PSF variability. We need to calibrate the PSF magnitudes separately from the + // Aperture-like magnitues. (We have an option to use the kron magnitudes or the + // other apertures here). I basically need to do this analysis separately for each + // magnitude type + + float MrelPSF = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); + if (isnan(MrelPSF)) { info->Nrel ++; continue; } + float MrelKron = getMrel (catalog, m, c, MAG_CLASS_KRON, MAG_SRC_CHP); // image.Mcal is not supposed to include the flat-field correction, so we need to @@ -1186,12 +1201,13 @@ // the flat-correction. NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat) - Mflat = getMflat (m, c, flatcorr, catalog); + float Mflat = getMflat (m, c, flatcorr, catalog); off_t n = catalog[c].measureT[m].averef; - Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); - if (isnan(Msys)) { + float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); + if (isnan(MsysPSF)) { info->Nsys++; continue; } + float MsysKron = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_KRON); PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode); @@ -1208,9 +1224,11 @@ assert (Nbright >= 0); - list[N] = Msys - Mrel - Mcal - Mgrid + Mflat; - dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); + psfMagList[N] = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat; + psfErrList[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); + kronMagList[N] = MsysKron - MrelKron - Mcal - Mgrid + Mflat; + kronErrList[N] = psfErrList[N]; if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) { - Mlist[Nbright] = list[N]; - dMlist[Nbright] = dlist[N]; + brightMagList[Nbright] = psfMagList[N]; + brightErrList[Nbright] = psfErrList[N]; Nbright ++; } @@ -1238,14 +1256,16 @@ } - liststats (list, dlist, NULL, N, &stats); + liststats (psfMagList, psfErrList, NULL, N, &stats); if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mmos: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N); // for now, I have no reason to measure these separately for camera-level images myMosaic[0].McalPSF = stats.mean; - myMosaic[0].McalAPER = stats.mean; myMosaic[0].dMcal = stats.error; myMosaic[0].McalChiSq = stats.chisq; myMosaic[0].nFitPhotom = N; + liststats (kronMagList, kronErrList, NULL, N, &stats); + myMosaic[0].McalAPER = stats.mean; + if (testImage) { fprintf (stderr, "test image %d (%d) %f %f %d ... ", (int) Nmos, myMosaic[0].start, stats.mean, stats.error, myMosaic[0].nFitPhotom); @@ -1254,9 +1274,9 @@ if (PLOTSTUFF) { fprintf (stderr, "Mmos: %6.3f %6.3f +/- %6.3f %5d %5d | %s\n", stats.mean, stats.median, stats.sigma, stats.Nmeas, N, image[MosaicToImage[Nmos][0]].name); - plot_setMcal (list, N, &stats, CLOUD_TOLERANCE); + plot_setMcal (psfMagList, N, &stats, CLOUD_TOLERANCE); } // bright end scatter - liststats (Mlist, dMlist, NULL, Nbright, &stats); + liststats (brightMagList, brightErrList, NULL, Nbright, &stats); myMosaic[0].dMsys = stats.sigma; @@ -1694,5 +1714,5 @@ c = MosaicToCatalog[i][j]; - Mcal = getMcal (m, c, flatcorr, catalog); + Mcal = getMcal (m, c, MAG_CLASS_PSF); if (isnan(Mcal)) continue; Mgrid = getMgrid (m, c); Index: /branches/eam_branches/ohana.20170822/src/relphot/src/StarOps.c =================================================================== --- /branches/eam_branches/ohana.20170822/src/relphot/src/StarOps.c (revision 40212) +++ /branches/eam_branches/ohana.20170822/src/relphot/src/StarOps.c (revision 40213) @@ -484,5 +484,5 @@ for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; - Mcal = getMcal (m, i, flatcorr, catalog); + Mcal = getMcal (m, i, MAG_CLASS_PSF); if (isnan(Mcal)) continue; Mmos = getMmos (m, i); @@ -688,5 +688,5 @@ // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers - Mcal = getMcal (m, i, flatcorr, catalog); + Mcal = getMcal (m, i, MAG_CLASS_PSF); if (isnan(Mcal)) { Ncal ++; continue; } Mmos = getMmos (m, i); @@ -739,5 +739,5 @@ // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers - Mcal = getMcal (m, i, flatcorr, catalog); + Mcal = getMcal (m, i, MAG_CLASS_PSF); if (isnan(Mcal)) continue; Mmos = getMmos (m, i); @@ -820,5 +820,5 @@ int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode); if (ecode != seccode) { continue;} - Mcal = getMcal (m, i, flatcorr, catalog); + Mcal = getMcal (m, i, MAG_CLASS_PSF); if (isnan(Mcal)) { continue;} Mmos = getMmos (m, i); Index: /branches/eam_branches/ohana.20170822/src/relphot/src/plot_scatter.c =================================================================== --- /branches/eam_branches/ohana.20170822/src/relphot/src/plot_scatter.c (revision 40212) +++ /branches/eam_branches/ohana.20170822/src/relphot/src/plot_scatter.c (revision 40213) @@ -40,5 +40,5 @@ if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; - Mcal = getMcal (m, i, flatcorr, catalog); + Mcal = getMcal (m, i, MAG_CLASS_PSF); if (isnan(Mcal)) continue; Mmos = getMmos (m, i); Index: /branches/eam_branches/ohana.20170822/src/relphot/src/setMrelCatalog.c =================================================================== --- /branches/eam_branches/ohana.20170822/src/relphot/src/setMrelCatalog.c (revision 40212) +++ /branches/eam_branches/ohana.20170822/src/relphot/src/setMrelCatalog.c (revision 40213) @@ -226,7 +226,6 @@ Mcal = measureT[k].Mcal; // check that this is zero for loaded REF value } else { - // getMcal_alt returns image[].Mcal modified by flatcorr(image.photom_map_id,x,y) - // NOTE: getMcal_alt does not include measure.Mflat - Mcal = getMcal_alt (meas, cat, flatcorr, measureT[k].Xccd, measureT[k].Yccd); + // getMcal returns image[].Mcal; note the flat-field correction is stored in measure.Mflat + Mcal = getMcal (meas, cat, MAG_CLASS_PSF); if (isnan(Mcal)) SKIP_THIS_MEAS(Ncal); Mmos = getMmos (meas, cat); @@ -530,5 +529,5 @@ off_t k; - float Mcal = 0, Mmos = 0, Mgrid = 0, Finst = 0; + float McalPSF = 0, McalAPER = 0, Mmos = 0, Mgrid = 0, Finst = 0; // set the primary projection cell and skycell for this coordinate @@ -630,8 +629,13 @@ // overlaps). Msys + measure.Mcal is our best guess of the true magnitude Mmos = Mgrid = 0; - Mcal = measure[k].Mcal; // check that this is zero for loaded REF value + McalPSF = measure[k].Mcal; // check that this is zero for loaded REF value + McalAPER = McalPSF; // check that this is zero for loaded REF value } else { - Mcal = getMcal_alt (meas, cat, flatcorr, measure[k].Xccd, measure[k].Yccd); - if (isnan(Mcal)) SKIP_THIS_MEAS_STACK(Ncal); + McalPSF = getMcal (meas, cat, MAG_CLASS_PSF); + if (isnan(McalPSF)) SKIP_THIS_MEAS_STACK(Ncal); // XXX really skip? + + McalAPER = getMcal (meas, cat, MAG_CLASS_PSF); + if (isnan(McalAPER)) SKIP_THIS_MEAS_STACK(Ncal); // XXX really skip? + Mmos = getMmos (meas, cat); if (isnan(Mmos)) SKIP_THIS_MEAS_STACK(Nmos); @@ -672,5 +676,7 @@ // get the zero point for the selected image - float zp = PhotZeroPoint (&measure[k], &average[0], &secfilt[0]) - (Mcal + Mmos + Mgrid); + // Use a different zero point for the PSF-like and APERTURE-like magnitudes + float zpPSF = PhotZeroPoint (&measure[k], &average[0], &secfilt[0]) - (McalPSF + Mmos + Mgrid); + float zpAPER = PhotZeroPoint (&measure[k], &average[0], &secfilt[0]) - (McalAPER + Mmos + Mgrid); // flux_cgs : erg sec^1 cm^-2 Hz^-1 @@ -704,17 +710,18 @@ // zpFactor to go from instrumental flux to Janskies - float zpFactor = pow(10.0, -0.4*zp + 3.56); + float zpFactorPSF = pow(10.0, -0.4*zpPSF + 3.56); + float zpFactorAPER = pow(10.0, -0.4*zpAPER + 3.56); // need to put in AB mag factor to get to Janskies (or uJy?) - secfilt[Nsec].FpsfStk = zpFactor * measure[k].FluxPSF; - secfilt[Nsec].dFpsfStk = zpFactor * measure[k].dFluxPSF; - secfilt[Nsec].FkronStk = zpFactor * measure[k].FluxKron; - secfilt[Nsec].dFkronStk = zpFactor * measure[k].dFluxKron; - secfilt[Nsec].FapStk = zpFactor * measure[k].FluxAp; + secfilt[Nsec].FpsfStk = zpFactorPSF * measure[k].FluxPSF; + secfilt[Nsec].dFpsfStk = zpFactorPSF * measure[k].dFluxPSF; + secfilt[Nsec].FkronStk = zpFactorAPER * measure[k].FluxKron; + secfilt[Nsec].dFkronStk = zpFactorAPER * measure[k].dFluxKron; + secfilt[Nsec].FapStk = zpFactorAPER * measure[k].FluxAp; // NOTE: for PV3, apFluxErr is broken (see pmSourcePhotometry.c:245). we are going to // use PSF flux error instead here: - // secfilt[Nsec].dFapStk = zpFactor * measure[k].dFluxAp; - secfilt[Nsec].dFapStk = zpFactor * measure[k].dFluxPSF; + // secfilt[Nsec].dFapStk = zpFactorAPER * measure[k].dFluxAp; + secfilt[Nsec].dFapStk = zpFactorAPER * measure[k].dFluxPSF; // Jy to AB mags @@ -859,7 +866,5 @@ Mcal = measure[k].Mcal; // check that this is zero for loaded REF value } else { - // use getMcal not getMcal_alt? - Mcal = getMcal_alt (meas, cat, NULL, measure[k].Xccd, measure[k].Yccd); - // Mcal = getMcal (meas, cat); + Mcal = getMcal (meas, cat, MAG_CLASS_PSF); if (isnan(Mcal)) SKIP_THIS_MEAS(Ncal); } Index: /branches/eam_branches/ohana.20170822/src/relphot/src/setMrelFinal.c =================================================================== --- /branches/eam_branches/ohana.20170822/src/relphot/src/setMrelFinal.c (revision 40212) +++ /branches/eam_branches/ohana.20170822/src/relphot/src/setMrelFinal.c (revision 40213) @@ -226,5 +226,5 @@ off_t Nim = getImageEntry (m, 0); if (Nim > -1) { - if (isnan(getMcal (m, 0, flatcorr, catalog))) goto skip; + if (isnan(getMcal (m, 0, MAG_CLASS_PSF))) goto skip; if (isnan(getMmos (m, 0))) goto skip; }