Index: trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- trunk/Ohana/src/relphot/include/relphot.h (revision 40019) +++ trunk/Ohana/src/relphot/include/relphot.h (revision 40291) @@ -77,11 +77,12 @@ unsigned int stop; short photcode; - float Mcal; + float McalPSF; + float McalAPER; float dMcal; float dMsys; - unsigned short nFitPhotom; - short Xm; + float McalChiSq; float secz; float ubercalDist; + unsigned short nFitPhotom; unsigned int flags; char skipCal; // if TRUE, this mosaic is incomplete and should not be calibrated @@ -201,8 +202,9 @@ typedef struct { - float Mcal; + float McalPSF; + float McalAPER; float dMcal; float dMagSys; - short Xm; + float McalChiSq; int nFitPhotom; int flags; @@ -220,5 +222,6 @@ typedef struct { - float Mcal; + float McalPSF; + float McalAPER; float dMcal; unsigned int imageID; @@ -336,4 +339,6 @@ int SyntheticPhotometry; +int USE_MCAL_PSF_FOR_STACK_APER; + char *PhotcodeList; @@ -404,4 +409,8 @@ int findMosaics PROTO((Catalog *catalog, int Ncatalog, int doMosaicList)); +void clearImages (void); +void checkImages (char *name); +int dumpMags (FILE *fout, Catalog *catalog, int Ncatalog); + void makeMosaics (Image *image, off_t Nimage, int mergeMcal); Mosaic *getMosaicForImage (off_t im); @@ -422,10 +431,9 @@ 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)); float getMmos PROTO((off_t meas, int cat)); -float getMrel PROTO((Catalog *catalog, off_t meas, int cat)); +float getMrel PROTO((Catalog *catalog, off_t meas, int cat, dvoMagClassType class, dvoMagSourceType source)); short getUbercalDist PROTO((off_t meas, int cat)); float getCenterOffset PROTO((off_t meas, int cat, Measure *measure, unsigned int *myID)); Index: trunk/Ohana/src/relphot/src/BrightCatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/BrightCatalog.c (revision 40019) +++ trunk/Ohana/src/relphot/src/BrightCatalog.c (revision 40291) @@ -66,5 +66,7 @@ GET_COLUMN(D, "DEC", double); GET_COLUMN(M, "MAG_SYS", float); - GET_COLUMN(Mcal, "MAG_CAL", float); + GET_COLUMN(Mkron, "MAG_KRON", float); + GET_COLUMN(McalPSF, "MCAL_PSF", float); + GET_COLUMN(McalAPER, "MCAL_APER", float); GET_COLUMN(Mflat, "MAG_FLAT", float); GET_COLUMN(dM, "MAG_ERR", float); @@ -89,5 +91,7 @@ measure[i].D = D[i]; measure[i].M = M[i]; - measure[i].Mcal = Mcal[i]; + measure[i].Mkron = Mkron[i]; + measure[i].McalPSF = McalPSF[i]; + measure[i].McalAPER = McalAPER[i]; measure[i].Mflat = Mflat[i]; measure[i].dM = dM[i]; @@ -110,5 +114,7 @@ free (D ); free (M ); - free (Mcal ); + free (Mkron ); + free (McalPSF ); + free (McalAPER); free (Mflat ); free (dM ); @@ -186,4 +192,5 @@ // need to create and assign to flat-field correction GET_COLUMN(M, "MAG", float); + GET_COLUMN(Mkron, "MAG_KRON", float); GET_COLUMN(dM, "MAG_ERR", float); GET_COLUMN(Mchisq, "MAG_CHI", float); @@ -199,16 +206,18 @@ ALLOCATE (secfilt, SecFilt, Nrow); for (i = 0; i < Nrow; i++) { - secfilt[i].M = M[i]; - secfilt[i].dM = dM[i]; - secfilt[i].Mchisq = Mchisq[i]; - secfilt[i].flags = flags[i]; - secfilt[i].Ncode = Ncode[i]; - secfilt[i].Nused = Nused[i]; - secfilt[i].Mmin = Mmin[i]; - secfilt[i].Mmax = Mmax[i]; + secfilt[i].MpsfChp = M[i]; + secfilt[i].MkronChp = Mkron[i]; + secfilt[i].dMpsfChp = dM[i]; + secfilt[i].Mchisq = Mchisq[i]; + secfilt[i].flags = flags[i]; + secfilt[i].Ncode = Ncode[i]; + secfilt[i].Nused = Nused[i]; + secfilt[i].Mmin = Mmin[i]; + secfilt[i].Mmax = Mmax[i]; } fprintf (stderr, "loaded data for %lld secfilt\n", (long long) Nrow); free (M ); + free (Mkron ); free (dM ); free (Mchisq); @@ -271,5 +280,7 @@ gfits_define_bintable_column (&theader, "D", "DEC", "dec", "degree", 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "MAG_SYS", "magnitude (sys)", NULL, 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "MAG_CAL", "magnitude (cal)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_KRON", "magnitude (sys,kron)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MCAL_PSF", "magnitude (cal)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MCAL_APER","magnitude (cal)", NULL, 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "MAG_FLAT", "magnitude (flat)", NULL, 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "MAG_ERR", "magnitude (err)", NULL, 1.0, 0.0); @@ -293,5 +304,7 @@ double *D ; ALLOCATE (D , double, catalog->Nmeasure); float *M ; ALLOCATE (M , float, catalog->Nmeasure); - float *Mcal ; ALLOCATE (Mcal , float, catalog->Nmeasure); + float *Mkron ; ALLOCATE (Mkron , float, catalog->Nmeasure); + float *McalPSF ; ALLOCATE (McalPSF , float, catalog->Nmeasure); + float *McalAPER ; ALLOCATE (McalAPER , float, catalog->Nmeasure); float *Mflat ; ALLOCATE (Mflat , float, catalog->Nmeasure); float *dM ; ALLOCATE (dM , float, catalog->Nmeasure); @@ -311,8 +324,10 @@ MeasureTiny *measure = catalog->measure; for (i = 0; i < catalog->Nmeasure; i++) { - R[i] = measure[i].R ; - D[i] = measure[i].D ; + R[i] = measure[i].R ; + D[i] = measure[i].D ; M[i] = measure[i].M ; - Mcal[i] = measure[i].Mcal ; + Mkron[i] = measure[i].Mkron ; + McalPSF[i] = measure[i].McalPSF ; + McalAPER[i] = measure[i].McalAPER ; Mflat[i] = measure[i].Mflat ; dM[i] = measure[i].dM ; @@ -334,5 +349,7 @@ gfits_set_bintable_column (&theader, &ftable, "DEC", D, catalog->Nmeasure); gfits_set_bintable_column (&theader, &ftable, "MAG_SYS", M, catalog->Nmeasure); - gfits_set_bintable_column (&theader, &ftable, "MAG_CAL", Mcal, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "MAG_KRON", Mkron, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "MCAL_PSF", McalPSF, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "MCAL_APER", McalAPER, catalog->Nmeasure); gfits_set_bintable_column (&theader, &ftable, "MAG_FLAT", Mflat, catalog->Nmeasure); gfits_set_bintable_column (&theader, &ftable, "MAG_ERR", dM, catalog->Nmeasure); @@ -352,5 +369,7 @@ free (D ); free (M ); - free (Mcal ); + free (Mkron ); + free (McalPSF ); + free (McalAPER); free (Mflat ); free (dM ); @@ -436,6 +455,7 @@ gfits_create_table_header (&theader, "BINTABLE", "SECFILT"); - gfits_define_bintable_column (&theader, "E", "MAG", "", "arcsec", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "MAG_ERR", "", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG", "", "mag", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_KRON", "", "mag", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_ERR", "", "mag", 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "MAG_CHI", "", NULL, 1.0, 0.0); gfits_define_bintable_column (&theader, "J", "FLAGS", "", NULL, 1.0, 0.0); @@ -454,4 +474,5 @@ // create intermediate storage arrays float *M ; ALLOCATE (M , float, Nsec); + float *Mkron ; ALLOCATE (Mkron , float, Nsec); float *dM ; ALLOCATE (dM , float, Nsec); float *Mchisq ; ALLOCATE (Mchisq , float, Nsec); @@ -465,16 +486,18 @@ SecFilt *secfilt = catalog->secfilt; for (i = 0; i < Nsec; i++) { - M [i] = secfilt[i].M ; - dM [i] = secfilt[i].dM ; - Mchisq[i] = secfilt[i].Mchisq ; - flags [i] = secfilt[i].flags ; - Ncode [i] = secfilt[i].Ncode ; - Nused [i] = secfilt[i].Nused ; - Mmin [i] = secfilt[i].Mmin ; - Mmax [i] = secfilt[i].Mmax ; + M [i] = secfilt[i].MpsfChp ; + Mkron [i] = secfilt[i].MkronChp; + dM [i] = secfilt[i].dMpsfChp; + Mchisq[i] = secfilt[i].Mchisq ; + flags [i] = secfilt[i].flags ; + Ncode [i] = secfilt[i].Ncode ; + Nused [i] = secfilt[i].Nused ; + Mmin [i] = secfilt[i].Mmin ; + Mmax [i] = secfilt[i].Mmax ; } // add the columns to the output array gfits_set_bintable_column (&theader, &ftable, "MAG", M , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_KRON", Mkron , Nsec); gfits_set_bintable_column (&theader, &ftable, "MAG_ERR", dM , Nsec); gfits_set_bintable_column (&theader, &ftable, "MAG_CHI", Mchisq, Nsec); @@ -486,4 +509,5 @@ free (M ); + free (Mkron ); free (dM ); free (Mchisq ); Index: trunk/Ohana/src/relphot/src/GridOps.c =================================================================== --- trunk/Ohana/src/relphot/src/GridOps.c (revision 40019) +++ trunk/Ohana/src/relphot/src/GridOps.c (revision 40291) @@ -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 ++; @@ -568,5 +568,5 @@ continue; } - Mrel = getMrel (catalog, m, c); + Mrel = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); if (isnan(Mrel)) { Nrel ++; @@ -645,9 +645,9 @@ continue; } - Mcal = getMcal (m, c, flatcorr, catalog); + Mcal = getMcal (m, c, MAG_CLASS_PSF); if (isnan(Mcal)) continue; Mmos = getMmos (m, c); if (isnan(Mmos)) continue; - Mrel = getMrel (catalog, m, c); + Mrel = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); if (isnan(Mrel)) continue; Index: trunk/Ohana/src/relphot/src/ImageMagIO.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageMagIO.c (revision 40019) +++ trunk/Ohana/src/relphot/src/ImageMagIO.c (revision 40291) @@ -58,12 +58,13 @@ char type[16]; - GET_COLUMN (Mcal, "MCAL", float); + GET_COLUMN (McalPSF, "MCAL_PSF", float); + GET_COLUMN (McalAPER, "MCAL_APER", float); GET_COLUMN (dMcal, "MCAL_ERR", float); GET_COLUMN (dMagSys, "MCAL_SYSERR", float); - GET_COLUMN (nFitPhotom, "NFIT", int); - GET_COLUMN (flags, "FLAGS", int); - GET_COLUMN (ubercalDist, "UDIST", int); - GET_COLUMN (imageID, "ID", int); - GET_COLUMN (Xm, "CHISQ", short); + GET_COLUMN (McalChiSq, "MCAL_CHISQ", float); + GET_COLUMN (nFitPhotom, "NFIT", int); + GET_COLUMN (flags, "FLAGS", int); + GET_COLUMN (ubercalDist, "UDIST", int); + GET_COLUMN (imageID, "ID", int); // free the memory associated with the FITS files @@ -75,8 +76,9 @@ ALLOCATE (image_mags, ImageMag, Nrow); for (i = 0; i < Nrow; i++) { - image_mags[i].Mcal = Mcal [i]; + image_mags[i].McalPSF = McalPSF [i]; + image_mags[i].McalAPER = McalAPER [i]; image_mags[i].dMcal = dMcal [i]; image_mags[i].dMagSys = dMagSys [i]; - image_mags[i].Xm = Xm [i]; + image_mags[i].McalChiSq = McalChiSq [i]; image_mags[i].nFitPhotom = nFitPhotom [i]; image_mags[i].flags = flags [i]; @@ -86,8 +88,9 @@ fprintf (stderr, "loaded data for %lld images\n", (long long) Nrow); - free (Mcal ); + free (McalPSF ); + free (McalAPER ); free (dMcal ); free (dMagSys ); - free (Xm ); + free (McalChiSq ); free (nFitPhotom ); free (flags ); @@ -121,20 +124,22 @@ gfits_create_table_header (&theader, "BINTABLE", "IMAGE_MAGS"); - gfits_define_bintable_column (&theader, "E", "MCAL", "cal offset", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MCAL_PSF", "PSF cal offset", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MCAL_APER", "APER cal offset", "magnitudes", 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "MCAL_ERR", "cal error", "magnitudes", 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "MCAL_SYSERR", "systematic error", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MCAL_CHISQ", "cal chisq", "unitless", 1.0, 0.0); gfits_define_bintable_column (&theader, "J", "NFIT", "number of fitted stars", "unitless", 1.0, 0.0); gfits_define_bintable_column (&theader, "J", "FLAGS", "analysis flags", "unitless", 1.0, 0.0); gfits_define_bintable_column (&theader, "J", "UDIST", "distance to ubercal images", "images", 1.0, 0.0); gfits_define_bintable_column (&theader, "J", "ID", "image ID", "unitless", 1.0, 0.0); - gfits_define_bintable_column (&theader, "I", "CHISQ", "cal chisq", "unitless", 1.0, FT_BZERO_INT16); // generate the output array that carries the data gfits_create_table (&theader, &ftable); - float *Mcal ; + float *McalPSF ; + float *McalAPER ; float *dMcal ; float *dMagSys ; - float *Xm ; + float *McalChiSq ; int *nFitPhotom ; int *flags ; @@ -143,8 +148,9 @@ // create intermediate storage arrays - ALLOCATE (Mcal , float, Nimage_mags); + ALLOCATE (McalPSF , float, Nimage_mags); + ALLOCATE (McalAPER , float, Nimage_mags); ALLOCATE (dMcal , float, Nimage_mags); ALLOCATE (dMagSys , float, Nimage_mags); - ALLOCATE (Xm , float, Nimage_mags); + ALLOCATE (McalChiSq , float, Nimage_mags); ALLOCATE (nFitPhotom , int, Nimage_mags); ALLOCATE (flags , int, Nimage_mags); @@ -154,8 +160,9 @@ // assign the storage arrays for (i = 0; i < Nimage_mags; i++) { - Mcal [i] = image_mags[i].Mcal ; + McalPSF [i] = image_mags[i].McalPSF ; + McalAPER [i] = image_mags[i].McalAPER ; dMcal [i] = image_mags[i].dMcal ; dMagSys [i] = image_mags[i].dMagSys ; - Xm [i] = image_mags[i].Xm ; + McalChiSq [i] = image_mags[i].McalChiSq ; nFitPhotom [i] = image_mags[i].nFitPhotom ; flags [i] = image_mags[i].flags ; @@ -165,8 +172,9 @@ // add the columns to the output array - gfits_set_bintable_column (&theader, &ftable, "MCAL", Mcal , Nimage_mags); + gfits_set_bintable_column (&theader, &ftable, "MCAL_PSF", McalPSF , Nimage_mags); + gfits_set_bintable_column (&theader, &ftable, "MCAL_APER", McalAPER , Nimage_mags); gfits_set_bintable_column (&theader, &ftable, "MCAL_ERR", dMcal , Nimage_mags); gfits_set_bintable_column (&theader, &ftable, "MCAL_SYSERR", dMagSys , Nimage_mags); - gfits_set_bintable_column (&theader, &ftable, "CHISQ", Xm , Nimage_mags); + gfits_set_bintable_column (&theader, &ftable, "MCAL_CHISQ", McalChiSq , Nimage_mags); gfits_set_bintable_column (&theader, &ftable, "NFIT", nFitPhotom , Nimage_mags); gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags , Nimage_mags); @@ -174,8 +182,9 @@ gfits_set_bintable_column (&theader, &ftable, "ID", imageID , Nimage_mags); - free (Mcal ); + free (McalPSF ); + free (McalAPER ); free (dMcal ); free (dMagSys ); - free (Xm ); + free (McalChiSq ); free (nFitPhotom ); free (flags ); Index: trunk/Ohana/src/relphot/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageOps.c (revision 40019) +++ trunk/Ohana/src/relphot/src/ImageOps.c (revision 40291) @@ -143,5 +143,6 @@ image[i].photom_map_id = input[i].photom_map_id; image[i].flags = input[i].flags ; - image[i].Mcal = input[i].Mcal ; + image[i].McalPSF = input[i].McalPSF ; + image[i].McalAPER = input[i].McalAPER ; image[i].dMcal = input[i].dMcal ; image[i].tzero = input[i].tzero ; @@ -179,5 +180,6 @@ subset[i].photom_map_id = image[i].photom_map_id; subset[i].flags = image[i].flags ; - subset[i].Mcal = image[i].Mcal ; + subset[i].McalPSF = image[i].McalPSF ; + subset[i].McalAPER = image[i].McalAPER ; subset[i].dMcal = image[i].dMcal ; subset[i].tzero = image[i].tzero ; @@ -430,58 +432,27 @@ } -// 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].Mcal; - - // 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].Mcal; - - // 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); -} - -// returns image.Mcal - ff(x,y) + 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 +} + short getUbercalDist (off_t meas, int cat) { @@ -499,5 +470,4 @@ } -// returns image.Mcal - ff(x,y) float getCenterOffset (off_t meas, int cat, Measure *measure, unsigned int *myID) { @@ -519,5 +489,4 @@ } -// returns image.Mcal - ff(x,y) int MatchImageName (off_t meas, int cat, char *name) { @@ -542,5 +511,4 @@ } -// returns image.Mcal - ff(x,y) int MatchImageSkycellID (off_t meas, int cat, int myTessID, int myProjectionID, int mySkycellID) { @@ -563,5 +531,4 @@ } -// returns image.Mcal - ff(x,y) int FindImageSkycellID (off_t meas, int cat, int *myTessID, int *myProjectionID, int *mySkycellID) { @@ -618,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; @@ -638,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); @@ -676,6 +643,7 @@ int minUbercalDist = 1000; - off_t Nref = 0; // total number of reference stars on the image - int Nbright = 0; // number of stars to measure the bright-end scatter + off_t Nref = 0; // number of stars used to measure McalPSF + int Nkron = 0; // number of stars to measure McalAPER + int Nbright = 0; // number of stars to measure the bright-end scatter if (N_onImage[i] == 0) { @@ -693,16 +661,23 @@ 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++; continue; } - Mrel = getMrel (catalog, m, c); - 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)) { Nrel ++; continue; @@ -714,9 +689,9 @@ // 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; @@ -731,19 +706,30 @@ 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"); } + float MrelKron = getMrel (catalog, m, c, MAG_CLASS_KRON, MAG_SRC_CHP); + float MsysKron = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_KRON); + + if (isfinite(MrelKron) && isfinite(MsysKron)) { + kronStars->flxlist[Nkron] = MsysKron - MrelKron - Mmos - Mgrid + Mflat; + kronStars->errlist[Nkron] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); + kronStars->wgtlist[Nkron] = 1; + kronStars->msklist[Nkron] = 0; + Nkron ++; + } + 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, MsysKron, MrelKron, Mmos, Mgrid, Mflat, kronStars->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; @@ -770,31 +756,39 @@ // 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].Mcal = stats.mean; + image[i].McalPSF = stats.mean; image[i].dMcal = stats.error; image[i].nFitPhotom = Nref; - image[i].Xm = 100.0*log10(stats.chisq); + image[i].McalChiSq = stats.chisq; 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, Nkron, &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, kronStars->flxlist[j], kronStars->errlist[j], kronStars->wgtlist[j], kronStars->msklist[j]); } } if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) { - fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].Mcal, image[i].dMcal); + fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].McalAPER, image[i].dMcal); } if (!mark && VERBOSE_IMAGE) { - fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].Mcal, image[i].dMcal); + fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].McalAPER, image[i].dMcal); } if (PLOTSTUFF) { - fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].Mcal, image[i].dMcal); - plot_setMcal (refStars->flxlist, Nref, &stats, CLOUD_TOLERANCE); + fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].McalAPER, image[i].dMcal); + plot_setMcal (psfStars->flxlist, Nref, &stats, CLOUD_TOLERANCE); } @@ -803,6 +797,6 @@ image[i].dMagSys = stats.sigma; - if (image[i].Mcal < -CLOUD_TOLERANCE) { - image[i].Mcal = 0.0; + if (image[i].McalPSF < -CLOUD_TOLERANCE) { + image[i].McalPSF = 0.0; } @@ -813,5 +807,6 @@ StatDataSetFree (brightStars, 1); - StatDataSetFree (refStars, 1); + StatDataSetFree (kronStars, 1); + StatDataSetFree (psfStars, 1); fprintf (stderr, "%d images calibrated\n", Ncalibrated); @@ -848,5 +843,5 @@ if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue; - mlist[N] = image[i].Mcal; + mlist[N] = image[i].McalPSF; slist[N] = image[i].dMcal; dlist[N] = 1; @@ -873,5 +868,5 @@ mark = FALSE; image[i].flags &= ~ID_IMAGE_PHOTOM_POOR; - mark = (image[i].dMcal > MaxScatter) || (fabs(image[i].Mcal - MedOffset) > MaxOffset); + mark = (image[i].dMcal > MaxScatter) || (fabs(image[i].McalPSF - MedOffset) > MaxOffset); if (mark) { Nmark ++; @@ -961,11 +956,11 @@ if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue; - Mlist[Nplot] = image[i].Mcal; + Mlist[Nplot] = image[i].McalPSF; dlist[Nplot] = image[i].dMcal; xlist[Nplot] = image[i].secz; minAirmass = MIN (image[i].secz, minAirmass); maxAirmass = MAX (image[i].secz, maxAirmass); - minMcal = MIN (image[i].Mcal, minMcal); - maxMcal = MAX (image[i].Mcal, maxMcal); + minMcal = MIN (image[i].McalPSF, minMcal); + maxMcal = MAX (image[i].McalPSF, maxMcal); mindMcal = MIN (image[i].dMcal, mindMcal); maxdMcal = MAX (image[i].dMcal, maxdMcal); @@ -1051,5 +1046,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); @@ -1092,5 +1087,5 @@ if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue; - list[n] = pow (10.0, 0.01*image[i].Xm); + list[n] = image[i].McalChiSq; dlist[n] = 1; n++; @@ -1125,5 +1120,5 @@ if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue; - list[n] = image[i].Mcal; + list[n] = image[i].McalPSF; dlist[n] = 1; n++; @@ -1170,2 +1165,21 @@ return (stats); } + +void clearImages (void) { + image = NULL; +} + +void checkImages (char *name) { + + off_t i; + fprintf (stderr, "--- %s ---\n", name); + + OhanaMemblock *ref = (OhanaMemblock *) image - 1; + fprintf (stderr, "file: %s\n", ref->file); + + fprintf (stderr, "0x%08llx \n", (unsigned long long) image); + + for (i = 0; i < Nimage; i++) { + fprintf (stderr, "%s : 0x%08x : %6.3f %6.3f\n", image[i].name, image[i].flags, image[i].McalPSF, image[i].McalAPER); + } +} Index: trunk/Ohana/src/relphot/src/ImageSubset.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageSubset.c (revision 40019) +++ trunk/Ohana/src/relphot/src/ImageSubset.c (revision 40291) @@ -63,5 +63,6 @@ char type[16]; - GET_COLUMN (Mcal, "MCAL", float); + GET_COLUMN (McalPSF, "MCAL_PSF", float); + GET_COLUMN (McalAPER, "MCAL_APER", float); GET_COLUMN (dMcal, "MCAL_ERR", float); GET_COLUMN (imageID, "IMAGE_ID", int); @@ -79,5 +80,6 @@ ALLOCATE (image, ImageSubset, Nrow); for (i = 0; i < Nrow; i++) { - image[i].Mcal = Mcal[i]; + image[i].McalPSF = McalPSF[i]; + image[i].McalAPER = McalAPER[i]; image[i].dMcal = dMcal[i]; image[i].imageID = imageID[i]; @@ -93,5 +95,6 @@ fprintf (stderr, "loaded data for %lld images\n", (long long) Nrow); - free (Mcal); + free (McalPSF); + free (McalAPER); free (dMcal); free (imageID); @@ -136,5 +139,6 @@ gfits_create_table_header (&theader, "BINTABLE", "IMAGE_SUBSET"); - gfits_define_bintable_column (&theader, "E", "MCAL", "zero point offset", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MCAL_PSF", "zero point offset", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MCAL_APER", "zero point offset", "magnitudes", 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "MCAL_ERR", "zero point error", "magnitudes", 1.0, 0.0); gfits_define_bintable_column (&theader, "J", "IMAGE_ID", "image ID", NULL, 1.0, FT_BZERO_INT32); @@ -151,5 +155,5 @@ gfits_create_table (&theader, &ftable); - float *Mcal, *dMcal; + float *McalPSF, *McalAPER, *dMcal; unsigned int *imageID, *map, *flags, *tzero; int *tessID, *projID, *skycellID; @@ -158,5 +162,6 @@ // create intermediate storage arrays - ALLOCATE (Mcal, float, Nimage); + ALLOCATE (McalPSF, float, Nimage); + ALLOCATE (McalAPER, float, Nimage); ALLOCATE (dMcal, float, Nimage); ALLOCATE (imageID, unsigned int, Nimage); @@ -172,5 +177,6 @@ // assign the storage arrays for (i = 0; i < Nimage; i++) { - Mcal[i] = image[i].Mcal; + McalPSF[i] = image[i].McalPSF; + McalAPER[i] = image[i].McalAPER; dMcal[i] = image[i].dMcal; imageID[i] = image[i].imageID; @@ -186,5 +192,6 @@ // add the columns to the output array - gfits_set_bintable_column (&theader, &ftable, "MCAL", Mcal, Nimage); + gfits_set_bintable_column (&theader, &ftable, "MCAL_PSF", McalPSF, Nimage); + gfits_set_bintable_column (&theader, &ftable, "MCAL_APER", McalAPER, Nimage); gfits_set_bintable_column (&theader, &ftable, "MCAL_ERR", dMcal, Nimage); gfits_set_bintable_column (&theader, &ftable, "IMAGE_ID", imageID, Nimage); @@ -198,5 +205,6 @@ gfits_set_bintable_column (&theader, &ftable, "UBERCAL_DIST", ucdist, Nimage); - free (Mcal); + free (McalPSF); + free (McalAPER); free (dMcal); free (imageID); Index: trunk/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- trunk/Ohana/src/relphot/src/MosaicOps.c (revision 40019) +++ trunk/Ohana/src/relphot/src/MosaicOps.c (revision 40291) @@ -111,13 +111,14 @@ /* a new mosaic, define ranges -- preserve the original values incase this image is not used */ - mosaic[Nmosaic].start = start; - mosaic[Nmosaic].stop = stop; - mosaic[Nmosaic].Mcal = 0.0; // note : mosaic stores only offsets relative to the original image values - mosaic[Nmosaic].dMcal = 0.0; // note : at the end, mosaic.Mcal is added back to the input images - mosaic[Nmosaic].dMsys = 0.0; - mosaic[Nmosaic].Xm = 0.0; - mosaic[Nmosaic].flags = image[i].flags; - mosaic[Nmosaic].secz = image[i].secz; - mosaic[Nmosaic].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode); + mosaic[Nmosaic].start = start; + mosaic[Nmosaic].stop = stop; + mosaic[Nmosaic].McalPSF = 0.0; // note : at the end, mosaic.Mcal is added back to the input images + mosaic[Nmosaic].McalAPER = 0.0; // note : mosaic stores only offsets relative to the original image values + mosaic[Nmosaic].dMcal = 0.0; // note : at the end, mosaic.Mcal is added back to the input images + mosaic[Nmosaic].dMsys = 0.0; + mosaic[Nmosaic].McalChiSq = 0.0;// NAN or 0.0? + mosaic[Nmosaic].flags = image[i].flags; + mosaic[Nmosaic].secz = image[i].secz; + mosaic[Nmosaic].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode); // XXX do we need to do something about flag consistency across a mosaic? @@ -208,14 +209,15 @@ for (i = 0; i < Nmosaic; i++) { /* a new mosaic, define ranges */ - mosaic[i].start = startTimesMosaic[i]; - mosaic[i].stop = 0; - mosaic[i].Mcal = 0.0; - mosaic[i].dMcal = 0.0; - mosaic[i].dMsys = 0.0; - mosaic[i].Xm = 0.0; - mosaic[i].flags = 0; - mosaic[i].secz = NAN; - mosaic[i].photcode = 0; - mosaic[i].skipCal = FALSE; + mosaic[i].start = startTimesMosaic[i]; + mosaic[i].stop = 0; + mosaic[i].McalPSF = 0.0; // note : at the end, mosaic.Mcal is added back to the input images + mosaic[i].McalAPER = 0.0; // note : mosaic stores only offsets relative to the original image values + mosaic[i].dMcal = 0.0; // note : at the end, mosaic.Mcal is added back to the input images + mosaic[i].dMsys = 0.0; + mosaic[i].McalChiSq = 0.0;// NAN or 0.0? + mosaic[i].flags = 0; + mosaic[i].secz = NAN; + mosaic[i].photcode = 0; + mosaic[i].skipCal = FALSE; memset (&mosaic[i].coords, 0, sizeof(Coords)); @@ -294,12 +296,15 @@ abort(); } - mosaic[j].stop = stop; - mosaic[j].Mcal = 0.0; - mosaic[j].dMcal = 0.0; - mosaic[j].Xm = 0.0; - mosaic[j].dMsys = subset[i].flags; - mosaic[j].flags = subset[i].flags; - mosaic[j].secz = subset[i].secz; - mosaic[j].photcode = GetPhotcodeEquivCodebyCode (subset[i].photcode); + mosaic[j].stop = stop; + + mosaic[j].McalPSF = 0.0; + mosaic[j].McalAPER = 0.0; + mosaic[j].dMcal = 0.0; + mosaic[j].McalChiSq = 0.0; + + mosaic[j].dMsys = subset[i].flags; + mosaic[j].flags = subset[i].flags; + mosaic[j].secz = subset[i].secz; + mosaic[j].photcode = GetPhotcodeEquivCodebyCode (subset[i].photcode); } @@ -377,14 +382,15 @@ // init the mosaic array values for (i = 0; i < Nmosaic; i++) { - mosaic[i].start = startTimesMosaic[i]; - mosaic[i].stop = 0; - mosaic[i].Mcal = 0.0; - mosaic[i].dMcal = 0.0; - mosaic[i].dMsys = 0.0; - mosaic[i].Xm = 0.0; - mosaic[i].flags = 0; - mosaic[i].secz = NAN; - mosaic[i].photcode = 0; - mosaic[i].skipCal = FALSE; + mosaic[i].start = startTimesMosaic[i]; + mosaic[i].stop = 0; + mosaic[i].McalPSF = 0.0; // note : at the end, mosaic.Mcal is added back to the input images + mosaic[i].McalAPER = 0.0; // note : mosaic stores only offsets relative to the original image values + mosaic[i].dMcal = 0.0; // note : at the end, mosaic.Mcal is added back to the input images + mosaic[i].dMsys = 0.0; + mosaic[i].McalChiSq = 0.0;// NAN or 0.0? + mosaic[i].flags = 0; + mosaic[i].secz = NAN; + mosaic[i].photcode = 0; + mosaic[i].skipCal = FALSE; memset (&mosaic[i].coords, 0, sizeof(Coords)); @@ -438,11 +444,12 @@ } mosaic[j].stop = stop; - mosaic[j].Mcal = 0.0; - mosaic[j].dMcal = 0.0; - mosaic[j].Xm = 0.0; - mosaic[j].dMsys = image[i].flags; - mosaic[j].flags = image[i].flags; - mosaic[j].secz = image[i].secz; - mosaic[j].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode); + mosaic[j].McalPSF = 0.0; + mosaic[j].McalAPER = 0.0; + mosaic[j].dMcal = 0.0; + mosaic[j].McalChiSq = 0.0; + mosaic[j].dMsys = image[i].flags; + mosaic[j].flags = image[i].flags; + mosaic[j].secz = image[i].secz; + mosaic[j].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode); } MARKTIME("assign images to mosaic: %f sec\n", dtime); @@ -594,7 +601,9 @@ mosaic[i].coords.cdelt2 = 1.0 / 3600.0; - mosaic[i].Mcal = 0.0; - mosaic[i].dMcal = 0.0; - mosaic[i].Xm = 0.0; + mosaic[i].McalPSF = 0.0; + mosaic[i].McalAPER = 0.0; + mosaic[i].dMcal = 0.0; + mosaic[i].dMsys = 0.0; + mosaic[i].McalChiSq = 0.0; } return; @@ -612,5 +621,5 @@ double dS, dX, dY; double R, D, Rmin, Rmax, Dmin, Dmax; - double Mcal, dMcal, Xm; + double McalPSF, McalAPER, dMcal, McalChiSq; fprintf (stderr, "*** moving Mcal from image.Mcal to mosaic.Mcal ***\n"); @@ -621,5 +630,5 @@ Dmax = Rmax = -360.0; dS = 0.0; - Mcal = dMcal = Xm = 0; + McalPSF = McalAPER = dMcal = McalChiSq = 0; for (j = 0; j < MosaicN_Image[i]; j++) { m = MosaicToImage[i][j]; @@ -678,11 +687,13 @@ // if it helps, note that ubercal uses a single zp per exposure, so the mean of those values is the same as the value - Mcal += image[m].Mcal; - dMcal += image[m].dMcal; - Xm += image[m].Xm; - - image[m].Mcal = 0.0; - image[m].dMcal = NAN; - image[m].Xm = NAN_S_SHORT; + McalPSF += image[m].McalPSF; + McalAPER += image[m].McalAPER; + dMcal += image[m].dMcal; + McalChiSq += image[m].McalChiSq; + + image[m].McalPSF = 0.0; + image[m].McalAPER = 0.0; + image[m].dMcal = NAN; + image[m].McalChiSq = NAN; } dS /= MosaicN_Image[i]; @@ -696,7 +707,8 @@ RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords); - mosaic[i].Mcal = Mcal / MosaicN_Image[i]; - mosaic[i].dMcal = dMcal / MosaicN_Image[i]; - mosaic[i].Xm = Xm / MosaicN_Image[i]; + mosaic[i].McalPSF = McalPSF / MosaicN_Image[i]; + mosaic[i].McalAPER = McalAPER / MosaicN_Image[i]; + mosaic[i].dMcal = dMcal / MosaicN_Image[i]; + mosaic[i].McalChiSq = McalChiSq / MosaicN_Image[i]; } if (!USE_GRID) return; @@ -724,14 +736,16 @@ for (j = 0; j < MosaicN_Image[i]; j++) { im = MosaicToImage[i][j]; - image[im].Mcal += mosaic[i].Mcal; - image[im].dMcal = mosaic[i].dMcal; - image[im].Xm = mosaic[i].Xm; + image[im].McalPSF += mosaic[i].McalPSF; + image[im].McalAPER += mosaic[i].McalAPER; + image[im].dMcal = mosaic[i].dMcal; + image[im].McalChiSq = mosaic[i].McalChiSq; image[im].ubercalDist = mosaic[i].ubercalDist; - image[im].dMagSys = mosaic[i].dMsys; - image[im].nFitPhotom = mosaic[i].nFitPhotom; - image[im].flags |= (mosaic[i].flags & ID_IMAGE_PHOTOM_FEW); - image[im].flags |= (mosaic[i].flags & ID_IMAGE_PHOTOM_POOR); - } - mosaic[i].Mcal = 0.0; + image[im].dMagSys = mosaic[i].dMsys; + image[im].nFitPhotom = mosaic[i].nFitPhotom; + image[im].flags |= (mosaic[i].flags & ID_IMAGE_PHOTOM_FEW); + image[im].flags |= (mosaic[i].flags & ID_IMAGE_PHOTOM_POOR); + } + mosaic[i].McalPSF = 0.0; + mosaic[i].McalAPER = 0.0; } } @@ -886,5 +900,5 @@ if (mosaic[i].flags & IMAGE_BAD) return (NAN); - value = mosaic[i].Mcal; + value = mosaic[i].McalPSF; return (value); } @@ -911,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; @@ -945,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); } @@ -1057,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); @@ -1125,23 +1147,23 @@ int N = 0; for (j = 0; j < N_onMosaic[Nmos]; j++) { - float Msys, Mrel, Mcal, Mgrid, Mflat; off_t m = MosaicToMeasure[Nmos][j]; off_t c = MosaicToCatalog[Nmos][j]; + // 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); - 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); } @@ -1150,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); - 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 @@ -1171,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); @@ -1193,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 ++; } @@ -1223,11 +1256,15 @@ } - 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); - myMosaic[0].Mcal = stats.mean; - myMosaic[0].dMcal = stats.error; + // for now, I have no reason to measure these separately for camera-level images + myMosaic[0].McalPSF = stats.mean; + myMosaic[0].dMcal = stats.error; + myMosaic[0].McalChiSq = stats.chisq; myMosaic[0].nFitPhotom = N; - myMosaic[0].Xm = 100.0*log10(stats.chisq); + + liststats (kronMagList, kronErrList, NULL, N, &stats); + myMosaic[0].McalAPER = stats.mean; if (testImage) { @@ -1237,17 +1274,17 @@ 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; - if (myMosaic[0].Mcal < -CLOUD_TOLERANCE) { - myMosaic[0].Mcal = 0.0; + if (myMosaic[0].McalPSF < -CLOUD_TOLERANCE) { + myMosaic[0].McalPSF = 0.0; } if (testImage) { - fprintf (stderr, "%f %f : %f\n", myMosaic[0].Mcal, myMosaic[0].dMsys, pow (10.0, 0.01*myMosaic[0].Xm)); + fprintf (stderr, "%f %f : %f\n", myMosaic[0].McalPSF, myMosaic[0].dMsys, myMosaic[0].McalChiSq); } @@ -1427,15 +1464,15 @@ imageOffset[i] = 0.0; - if (VERBOSE2 && (fabs(mosaic[i].Mcal) < CLOUD_TOLERANCE)) { - fprintf (stderr, "cloud-free: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].Mcal); + if (VERBOSE2 && (fabs(mosaic[i].McalPSF) < CLOUD_TOLERANCE)) { + fprintf (stderr, "cloud-free: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].McalPSF); } - if (VERBOSE2 && (mosaic[i].Mcal < -CLOUD_TOLERANCE)) { - imageOffset[i] = -mosaic[i].Mcal; - // NOTE the negative sign: down below, we are going to add in the negative of Mcal + if (VERBOSE2 && (mosaic[i].McalPSF < -CLOUD_TOLERANCE)) { + imageOffset[i] = -mosaic[i].McalPSF; + // NOTE the negative sign: down below, we are going to add in the negative of McalPSF // to this image, and the propagated mean values for other images - fprintf (stderr, "anti-clouds: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].Mcal); + fprintf (stderr, "anti-clouds: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].McalPSF); } - if (VERBOSE2 && (mosaic[i].Mcal > CLOUD_TOLERANCE)) { - fprintf (stderr, "cloudy : %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].Mcal); + if (VERBOSE2 && (mosaic[i].McalPSF > CLOUD_TOLERANCE)) { + fprintf (stderr, "cloudy : %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].McalPSF); } } @@ -1540,5 +1577,5 @@ } - // a negative cloud image (cloud: Mcal > 0; anti-clouds: Mcal < 0; imageOffset = -Mcal) + // a negative cloud image (cloud: McalPSF > 0; anti-clouds: McalPSF < 0; imageOffset = -McalPSF) if (imageOffset[i] > 0.0) continue; @@ -1564,5 +1601,5 @@ // find the images / mosaics with negative clouds and save their offset for (i = 0; i < Nmosaic; i++) { - mosaic[i].Mcal += imageOffset[i]; + mosaic[i].McalPSF += imageOffset[i]; } @@ -1605,5 +1642,5 @@ if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; - list[n] = mosaic[i].Mcal; + list[n] = mosaic[i].McalPSF; dlist[n] = 1; n++; @@ -1677,9 +1714,9 @@ 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); if (isnan(Mgrid)) continue; - Mrel = getMrel (catalog, m, c); + Mrel = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); if (isnan(Mrel)) continue; N++; @@ -1719,5 +1756,5 @@ if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; - list[n] = pow (10.0, 0.01*mosaic[i].Xm); + list[n] = mosaic[i].McalChiSq; dlist[n] = 1; n++; @@ -1751,5 +1788,5 @@ if (mosaic[i].flags & IMAGE_BAD) continue; if (mosaic[i].skipCal) continue; - mlist[N] = mosaic[i].Mcal; + mlist[N] = mosaic[i].McalPSF; slist[N] = mosaic[i].dMcal; dlist[N] = 1; @@ -1781,5 +1818,5 @@ Nscatter ++; } - if (fabs(mosaic[i].Mcal - MedOffset) > MaxOffset) { + if (fabs(mosaic[i].McalPSF - MedOffset) > MaxOffset) { mark = TRUE; Noffset ++; @@ -1858,5 +1895,5 @@ for (i = 0; i < Nmosaic; i++) { - Mlist[i] = mosaic[i].Mcal; + Mlist[i] = mosaic[i].McalPSF; dlist[i] = mosaic[i].dMcal; xlist[i] = mosaic[i].secz; Index: trunk/Ohana/src/relphot/src/StarOps.c =================================================================== --- trunk/Ohana/src/relphot/src/StarOps.c (revision 40019) +++ trunk/Ohana/src/relphot/src/StarOps.c (revision 40291) @@ -34,22 +34,66 @@ } -float getMrel (Catalog *catalog, off_t meas, int cat) { - - int Nsec, Nsecfilt, photcode; - int ave; - float value; - - ave = catalog[cat].measureT[meas].averef; - photcode = catalog[cat].measureT[meas].photcode; +// +float getMrel (Catalog *catalog, off_t meas, int cat, dvoMagClassType class, dvoMagSourceType source) { + + int ave = catalog[cat].measureT[meas].averef; + int photcode = catalog[cat].measureT[meas].photcode; int ecode = GetPhotcodeEquivCodebyCode (photcode); - Nsec = GetPhotcodeNsec(ecode); - Nsecfilt = GetPhotcodeNsecfilt (); + int Nsec = GetPhotcodeNsec(ecode); + int Nsecfilt = GetPhotcodeNsecfilt (); + + int entry = Nsecfilt*ave+Nsec; + + SecFilt *secfilt = &catalog[cat].secfilt[entry]; // is this star OK? - if (catalog[cat].secfilt[Nsecfilt*ave+Nsec].flags & STAR_BAD) return (NAN); - - value = catalog[cat].secfilt[Nsecfilt*ave+Nsec].M; - return (value); + if (secfilt->flags & STAR_BAD) return (NAN); + + switch (class) { + case MAG_CLASS_PSF: + switch (source) { + case MAG_SRC_CHP: + return secfilt->MpsfChp; + case MAG_SRC_WRP: + return secfilt->MpsfWrp; + case MAG_SRC_STK: + return secfilt->MpsfStk; + default: + return NAN; + } + break; + case MAG_CLASS_KRON: + switch (source) { + case MAG_SRC_CHP: + return secfilt->MkronChp; + case MAG_SRC_WRP: + return secfilt->MkronWrp; + case MAG_SRC_STK: + return secfilt->MkronStk; + default: + return NAN; + } + break; + case MAG_CLASS_APER: + switch (source) { + case MAG_SRC_CHP: + return secfilt->MapChp; + case MAG_SRC_WRP: + return secfilt->MapWrp; + case MAG_SRC_STK: + return secfilt->MapStk; + default: + return NAN; + } + break; + case MAG_CLASS_NONE: + case MAG_CLASS_DEV: // DeVaucouleur Model (only for galphot) + case MAG_CLASS_EXP: // Exponential Model (only for galphot) + case MAG_CLASS_SER: // Sersic Model (only for galphot) + return NAN; + break; + } + return NAN; // should not be able to reach here } @@ -429,5 +473,4 @@ int i; off_t j, k, m; - float Mcal, Mmos, Mgrid; MEAS_BAD = ID_MEAS_NOCAL; @@ -440,10 +483,18 @@ 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); - if (isnan(Mcal)) continue; - Mmos = getMmos (m, i); + float McalPSF = getMcal (m, i, MAG_CLASS_PSF); + float McalAPER = getMcal (m, i, MAG_CLASS_APER); + if (isnan(McalPSF)) continue; + float Mmos = getMmos (m, i); if (isnan(Mmos)) continue; - Mgrid = getMgrid (m, i); + float Mgrid = getMgrid (m, i); if (isnan(Mgrid)) continue; + + // Note that this operation is setting measure->McalAPER to image->McalAPER only + // for the stacks. Elsewhere (setMrelCatalog) we are using image->McalPSF for all + // types of measurements EXCEPT stacks. + + // we need to use McalAPER for stacks (and only stacks) + int useStackAper = !USE_MCAL_PSF_FOR_STACK_APER && isGPC1stack(catalog[i].measureT[m].photcode); // note that measurements for which the image is not selected will not be modified @@ -451,8 +502,9 @@ // set the output calibration - catalog[i].measure[m].Mcal = Mcal + Mmos + Mgrid; + catalog[i].measure[m].McalPSF = McalPSF + Mmos + Mgrid; + catalog[i].measure[m].McalAPER = useStackAper ? McalAPER + Mmos + Mgrid : McalPSF + Mmos + Mgrid; if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) { - myAssert (isfinite(catalog[i].measure[m].Mcal), "oops, broke an ubercal mag"); + myAssert (isfinite(catalog[i].measure[m].McalPSF), "oops, broke an ubercal mag"); } } @@ -501,4 +553,26 @@ } +int dumpMags (FILE *fout, Catalog *catalog, int Ncatalog) { + + int i, n; + off_t j; + + int Nsecfilt = GetPhotcodeNsecfilt (); + + for (i = 0; i < Ncatalog; i++) { + for (j = 0; (j < catalog[i].Naverage) && (j < 2); j++) { + fprintf (fout, "%08x %08x %10.6f %10.6f : \n", catalog[i].averageT[j].catID, catalog[i].averageT[j].objID, catalog[i].averageT[j].R, catalog[i].averageT[j].D); + for (n = 0; n < 5; n++) { + fprintf (fout, "secf %5d | %6.3f %6.3f | %6.3f %6.3f | %6.3f %6.3f\n", n, catalog[i].secfilt[j*Nsecfilt + n].MpsfChp, catalog[i].secfilt[j*Nsecfilt + n].MkronChp, catalog[i].secfilt[j*Nsecfilt + n].MpsfStk, catalog[i].secfilt[j*Nsecfilt + n].MkronStk, catalog[i].secfilt[j*Nsecfilt + n].MpsfWrp, catalog[i].secfilt[j*Nsecfilt + n].MkronWrp); + } + int m = catalog[i].averageT[j].measureOffset; + for (n = 0; n < catalog[i].averageT[j].Nmeasure; n++) { + fprintf (fout, "meas %5d %5d | %6.3f %6.3f | %6.3f %6.3f\n", m+n, catalog[i].measureT[m+n].photcode, catalog[i].measureT[m+n].M, catalog[i].measureT[m+n].Mkron, catalog[i].measureT[m+n].McalPSF, catalog[i].measureT[m+n].McalAPER); + } + } + } + return TRUE; +} + void clean_stars (Catalog *catalog, int Ncatalog) { @@ -513,5 +587,5 @@ if (VERBOSE) fprintf (stderr, "marking poor stars\n"); - /* find Xm median -> ChiSq lim must be > median */ + /* find Mchisq median -> ChiSq lim must be > median */ for (i = Ntot = 0; i < Ncatalog; i++) { Ntot += catalog[i].Naverage; @@ -535,5 +609,5 @@ if (isnan(Mchisq)) continue; xlist[Ntot] = Mchisq; - slist[Ntot] = catalog[i].secfilt[Nsecfilt*j+Nsec].dM; + slist[Ntot] = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp; dlist[Ntot] = 1; Ntot ++; @@ -551,5 +625,5 @@ for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Naverage; j++) { - dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dM; + dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp; float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq; mark = (dM > MaxScatter) || (isnan(Mchisq)) || (Mchisq > MaxChisq); @@ -578,4 +652,6 @@ // NSIGMA_REJECT (5) sigma of the mean +// this function only operations on the PSF magnitudes + # define NSIGMA_CLIP 3.0 # define NSIGMA_REJECT 5.0 @@ -644,5 +720,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); @@ -695,5 +771,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); @@ -776,5 +852,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); @@ -865,5 +941,5 @@ if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; - dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dM; + dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp; if (isnan(dM)) continue; list[n] = dM; @@ -906,5 +982,5 @@ for (j = 0; j < catalog[i].Naverage; j++) { if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; - dMrel = catalog[i].secfilt[Nsecfilt*j+Nsec].dM; + dMrel = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp; bin = dMrel / 0.00025; bin = MAX (0, MIN (NBIN-1, bin)); @@ -946,5 +1022,5 @@ for (j = 0; j < catalog[i].Naverage; j++) { if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; - xlist[N] = catalog[i].secfilt[Nsecfilt*j+Nsec].M; + xlist[N] = catalog[i].secfilt[Nsecfilt*j+Nsec].MpsfChp; value = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq; if (isnan((double)(value))) continue; Index: trunk/Ohana/src/relphot/src/args.c =================================================================== --- trunk/Ohana/src/relphot/src/args.c (revision 40019) +++ trunk/Ohana/src/relphot/src/args.c (revision 40291) @@ -325,4 +325,9 @@ remove_argument (N, &argc, argv); CALIBRATE_STACKS_AND_WARPS = TRUE; + } + USE_MCAL_PSF_FOR_STACK_APER = FALSE; + if ((N = get_argument (argc, argv, "-use-mcal-psf-for-stack-aper"))) { + remove_argument (N, &argc, argv); + USE_MCAL_PSF_FOR_STACK_APER = TRUE; } @@ -696,4 +701,10 @@ } + USE_MCAL_PSF_FOR_STACK_APER = FALSE; + if ((N = get_argument (argc, argv, "-use-mcal-psf-for-stack-aper"))) { + remove_argument (N, &argc, argv); + USE_MCAL_PSF_FOR_STACK_APER = TRUE; + } + /* define time */ TimeSelect = FALSE; Index: trunk/Ohana/src/relphot/src/bcatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/bcatalog.c (revision 40019) +++ trunk/Ohana/src/relphot/src/bcatalog.c (revision 40291) @@ -153,10 +153,12 @@ if (subcatalog[0].measureT[Nmeasure].dbFlags & ID_MEAS_PHOTOM_UBERCAL) { if (!KEEP_UBERCAL) { - subcatalog[0].measureT[Nmeasure].Mcal = 0.0; + subcatalog[0].measureT[Nmeasure].McalPSF = 0.0; + subcatalog[0].measureT[Nmeasure].McalAPER = 0.0; subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL; } } else { if (getImageEntry (Nmeasure, Ncat) >= 0) { - subcatalog[0].measureT[Nmeasure].Mcal = 0.0; + subcatalog[0].measureT[Nmeasure].McalPSF = 0.0; + subcatalog[0].measureT[Nmeasure].McalAPER = 0.0; } } Index: trunk/Ohana/src/relphot/src/launch_region_hosts.c =================================================================== --- trunk/Ohana/src/relphot/src/launch_region_hosts.c (revision 40019) +++ trunk/Ohana/src/relphot/src/launch_region_hosts.c (revision 40291) @@ -58,22 +58,25 @@ strextend (&command, "-cloud-limit %f", CLOUD_TOLERANCE); - if (VERBOSE) strextend (&command, "-v"); - if (VERBOSE2) strextend (&command, "-vv"); - if (RESET) strextend (&command, "-reset"); - if (RESET_ZEROPTS) strextend (&command, "-reset-zpts"); - if (!KEEP_UBERCAL) strextend (&command, "-reset-ubercal"); - if (DophotSelect) strextend (&command, "-dophot %d", DophotValue); - if (ImagSelect) strextend (&command, "-instmag %f %f", ImagMin, ImagMax); - if (MaxDensityUse) strextend (&command, "-max-density %f", MaxDensityValue); - if (SyntheticPhotometry) strextend (&command, "-synthphot"); - if (USE_BASIC_CHECK) strextend (&command, "-basic-image-search"); + if (VERBOSE) strextend (&command, "-v"); + if (VERBOSE2) strextend (&command, "-vv"); + if (RESET) strextend (&command, "-reset"); + if (RESET_ZEROPTS) strextend (&command, "-reset-zpts"); + if (!KEEP_UBERCAL) strextend (&command, "-reset-ubercal"); + if (DophotSelect) strextend (&command, "-dophot %d", DophotValue); + if (ImagSelect) strextend (&command, "-instmag %f %f", ImagMin, ImagMax); + if (MaxDensityUse) strextend (&command, "-max-density %f", MaxDensityValue); + if (SyntheticPhotometry) strextend (&command, "-synthphot"); + if (USE_BASIC_CHECK) strextend (&command, "-basic-image-search"); - if (UPDATE) strextend (&command, "-update"); - if (MOSAIC_ZEROPT) strextend (&command, "-mosaic"); - if (FREEZE_IMAGES) strextend (&command, "-imfreeze"); - if (FREEZE_MOSAICS) strextend (&command, "-mosfreeze"); - if (PARALLEL) strextend (&command, "-parallel"); - if (PARALLEL_MANUAL) strextend (&command, "-parallel-manual"); - if (PARALLEL_SERIAL) strextend (&command, "-parallel-serial"); + if (UPDATE) strextend (&command, "-update"); + if (MOSAIC_ZEROPT) strextend (&command, "-mosaic"); + if (FREEZE_IMAGES) strextend (&command, "-imfreeze"); + if (FREEZE_MOSAICS) strextend (&command, "-mosfreeze"); + if (CALIBRATE_STACKS_AND_WARPS) strextend (&command, "-only-stacks-and-warps"); + if (USE_MCAL_PSF_FOR_STACK_APER) { strextend (&command, "-use-mcal-psf-for-stack-aper"); } + + if (PARALLEL) strextend (&command, "-parallel"); + if (PARALLEL_MANUAL) strextend (&command, "-parallel-manual"); + if (PARALLEL_SERIAL) strextend (&command, "-parallel-serial"); // XXX deprecate this if we are happy with the new version Index: trunk/Ohana/src/relphot/src/liststats.c =================================================================== --- trunk/Ohana/src/relphot/src/liststats.c (revision 40019) +++ trunk/Ohana/src/relphot/src/liststats.c (revision 40291) @@ -194,5 +194,5 @@ # define WEIGHT_THRESHOLD 0.3 -int fit_least_squares (double *fit, double *y, double *dy, double *wgt, double *wt, int Npts); +int fit_least_squares (double *fit, double *err, double *y, double *dy, double *wgt, double *wt, int Npts); double VectorFractionInterpolate (double *values, float fraction, int Npts); double weight_cauchy (double x); @@ -253,5 +253,5 @@ double oldValue = value; - if (!fit_least_squares (&value, dataset->flxlist, dataset->errlist, dataset->wgtlist, dataset->wtvals, Npoints)) { + if (!fit_least_squares (&value, NULL, dataset->flxlist, dataset->errlist, dataset->wgtlist, dataset->wtvals, Npoints)) { value = oldValue; break; @@ -330,5 +330,5 @@ } - if (!fit_least_squares (&value, dataset->ysample, dataset->dysample, dataset->wtsample, NULL, Nkeep)) continue; + if (!fit_least_squares (&value, NULL, dataset->ysample, dataset->dysample, dataset->wtsample, NULL, Nkeep)) continue; dataset->bvalue[Nboot] = value; @@ -342,9 +342,16 @@ stats->error = (Shi - Slo) / 2.0; + // bootstrap can sometimes yield an excessively-optimistic result for the error. Do not let + // the reported error be smaller than the formal error + double errvalue; + if (fit_least_squares (&value, &errvalue, dataset->ykeep, dataset->dykeep, dataset->wtkeep, NULL, Nkeep)) { + stats->error = MAX (stats->error, errvalue); + } + return TRUE; } // wgt is externally-supplied weight, wt is optional -int fit_least_squares (double *fit, double *y, double *dy, double *wgt, double *wt, int Npts) { +int fit_least_squares (double *fit, double *err, double *y, double *dy, double *wgt, double *wt, int Npts) { int i; @@ -364,5 +371,6 @@ } if (S0 == 0.0) return FALSE; - *fit = S1 / S0; + *fit = S1 / S0; + if (err) { *err = 1.0 / S0; } return TRUE; } Index: trunk/Ohana/src/relphot/src/load_images.c =================================================================== --- trunk/Ohana/src/relphot/src/load_images.c (revision 40019) +++ trunk/Ohana/src/relphot/src/load_images.c (revision 40291) @@ -52,4 +52,24 @@ free (subset); subset = image; + } + } + + // reset image values as needed. always allow 'few' images to succeed, if possible (new + // images / detections may have been added + for (off_t i = 0; i < Nsubset; i++) { + if (subset[i].flags & ID_IMAGE_PHOTOM_FEW) { + subset[i].flags &= ~(ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR); + } + if (RESET) { + if (RESET_ZEROPTS) { + if (!KEEP_UBERCAL || !(subset[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) { + subset[i].McalPSF = 0.0; + subset[i].McalAPER = 0.0; + subset[i].dMcal = NAN; + subset[i].flags &= ~ID_IMAGE_PHOTOM_UBERCAL; + } + } + subset[i].flags &= ~ID_IMAGE_PHOTOM_POOR; + subset[i].ubercalDist = 1000; } } Index: trunk/Ohana/src/relphot/src/plot_scatter.c =================================================================== --- trunk/Ohana/src/relphot/src/plot_scatter.c (revision 40019) +++ trunk/Ohana/src/relphot/src/plot_scatter.c (revision 40291) @@ -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); @@ -47,5 +47,5 @@ if (isnan(Mgrid)) continue; - Mrel = catalog[i].secfilt[Nsecfilt*j+Nsec].M; + Mrel = catalog[i].secfilt[Nsecfilt*j+Nsec].MpsfChp; if (isnan(Mrel)) continue; Index: trunk/Ohana/src/relphot/src/plotstuff.c =================================================================== --- trunk/Ohana/src/relphot/src/plotstuff.c (revision 40019) +++ trunk/Ohana/src/relphot/src/plotstuff.c (revision 40291) @@ -179,5 +179,5 @@ KapaInitGraph (graphdata); - graphdata[0].style = 2; + graphdata[0].style = KAPA_PLOT_POINTS; /* points */ graphdata[0].ptype = 2; graphdata[0].ltype = 0; Index: trunk/Ohana/src/relphot/src/reload_catalogs.c =================================================================== --- trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 40019) +++ trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 40291) @@ -247,4 +247,5 @@ if (SYNTH_ZERO_POINTS) { strextend (&command, "-synthphot-zpts %s", SYNTH_ZERO_POINTS); } if (USE_BASIC_CHECK) { strextend (&command, "-basic-image-search"); } + if (USE_MCAL_PSF_FOR_STACK_APER) { strextend (&command, "-use-mcal-psf-for-stack-aper"); } if (!(STAGES & STAGE_CHIP)) { strextend (&command, "-skip-chip"); } Index: trunk/Ohana/src/relphot/src/relphot_images.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_images.c (revision 40019) +++ trunk/Ohana/src/relphot/src/relphot_images.c (revision 40291) @@ -42,8 +42,11 @@ MARKTIME("-- load images: %f sec\n", dtime); + // checkImages("load images"); + /* unlock, if we can (else, unlocked below) */ if (!UPDATE) dvo_image_unlock (&db); // load the flat correction table (if defined) + // we no longer need to do this: the flatcorrection is now saved in measure.Mflat # if (0) char flatcorrfile[256]; @@ -60,4 +63,6 @@ /* add in a loop over the catalogs calling dvo_catalog_chipcoords */ + // checkImages("load catalog"); + /* match measurements with images, mosaics */ initImageBins (catalog, Ncatalog, TRUE); @@ -77,4 +82,5 @@ setExclusions (catalog, Ncatalog, TRUE); + // checkImages("set Exclu"); global_stats (catalog, Ncatalog, flatcorr, 0); @@ -108,9 +114,18 @@ plot_scatter (catalog, Ncatalog, flatcorr); } + if (VERBOSE2) dumpMags (stderr, catalog, Ncatalog); + // checkImages("start loop"); + setMrel (catalog, Ncatalog, flatcorr); // threaded + if (VERBOSE2) dumpMags (stderr, catalog, Ncatalog); + // checkImages("set Mrel"); + if (PLOTSTUFF) { plot_scatter (catalog, Ncatalog, flatcorr); } setMcal (catalog, FALSE, flatcorr); + // checkImages("set Mcal"); + if (VERBOSE2) dumpMags (stderr, catalog, Ncatalog); + setMmos (catalog, FALSE, flatcorr); setMgrid (catalog, flatcorr); @@ -213,4 +228,6 @@ // ensure the db format is updated dbX.format = dvo_catalog_catformat (UPDATE_CATFORMAT); + gfits_modify (&dbX.header, "FORMAT", "%s", 1, UPDATE_CATFORMAT); + char photcodeFile[1024]; sprintf (photcodeFile, "%s/Photcodes.dat", CATDIR); @@ -257,4 +274,5 @@ // ensure the db format is updated db.format = dvo_catalog_catformat (UPDATE_CATFORMAT); + gfits_modify (&db.header, "FORMAT", "%s", 1, UPDATE_CATFORMAT); } if (CALIBRATE_STACKS_AND_WARPS || (NLOOP > 0)) { Index: trunk/Ohana/src/relphot/src/relphot_objects.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_objects.c (revision 40019) +++ trunk/Ohana/src/relphot/src/relphot_objects.c (revision 40291) @@ -251,4 +251,5 @@ if (USE_BASIC_CHECK) { strextend (&command, "-basic-image-search"); } if (USE_ALL_IMAGES) { strextend (&command, "-use-all-images"); } + if (USE_MCAL_PSF_FOR_STACK_APER) { strextend (&command, "-use-mcal-psf-for-stack-aper"); } if (!(STAGES & STAGE_CHIP)) { strextend (&command, "-skip-chip"); } Index: trunk/Ohana/src/relphot/src/relphot_parallel_images.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_parallel_images.c (revision 40019) +++ trunk/Ohana/src/relphot/src/relphot_parallel_images.c (revision 40291) @@ -36,7 +36,4 @@ client_logger_message ("loaded images\n"); - // once we have read this table, we should remove it for repeat runs - // unlink (IMAGE_TABLE); // XXX a bit risky, add some protection? - makeMosaics (image, Nimage, TRUE); @@ -44,4 +41,5 @@ // load the flat correction table (if defined) + // we no longer need to do this: the flatcorrection is now saved in measure.Mflat # if (0) char flatcorrfile[256]; @@ -150,7 +148,8 @@ client_logger_message ("all hosts are done the loops\n"); - /* set Mcal & Mmos for bad images */ - setMcal (catalog, TRUE, flatcorr); - setMmos (catalog, TRUE, flatcorr); + /* set Mcal & Mmos for bad images (for stack_and_warps, force images to be measured) */ + int onlyPoorImages = !CALIBRATE_STACKS_AND_WARPS; + setMcal (catalog, onlyPoorImages, flatcorr); + setMmos (catalog, onlyPoorImages, flatcorr); MARKTIME("-- finalize Mcal values: %f sec\n", dtime); @@ -171,4 +170,5 @@ freeGridBins (Ncatalog); freeImages((char *)image); + free (image); SkyListFree(skylist); Index: trunk/Ohana/src/relphot/src/relphot_parallel_regions.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_parallel_regions.c (revision 40019) +++ trunk/Ohana/src/relphot/src/relphot_parallel_regions.c (revision 40291) @@ -51,4 +51,5 @@ // ensure the db format is updated db.format = dvo_catalog_catformat (UPDATE_CATFORMAT); + gfits_modify (&db.header, "FORMAT", "%s", 1, UPDATE_CATFORMAT); } @@ -84,5 +85,9 @@ // save the changes to the image parameters - dvo_image_save (&db, VERBOSE); // this function modifies the db.ftable.buffer: do not free stored Image table + dvo_image_save (&db, VERBOSE); + + // dvo_image_save frees db.ftable.buffer (== image) and replaces it: do not free stored Image table + clearImages (); + dvo_image_unlock (&db); MARKTIME ("finished relphot -parallel-regions: %f sec total\n", dtime); Index: trunk/Ohana/src/relphot/src/relphot_synthphot_catalog.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_synthphot_catalog.c (revision 40019) +++ trunk/Ohana/src/relphot/src/relphot_synthphot_catalog.c (revision 40291) @@ -92,8 +92,8 @@ float *value = (float *) zpts->matrix[Nsec].buffer; float ZP = !isnan(value[Npix]) ? value[Npix] : 0.0; - secfilt[Nsec].M = measure[measSYN[i]].M + ZP; - secfilt[Nsec].dM = 0.6; - secfilt[Nsec].Mchisq = NAN; - secfilt[Nsec].flags |= ID_SECF_USE_SYNTH; + secfilt[Nsec].MpsfChp = measure[measSYN[i]].M + ZP; + secfilt[Nsec].dMpsfChp = 0.6; + secfilt[Nsec].Mchisq = NAN; + secfilt[Nsec].flags |= ID_SECF_USE_SYNTH; } } Index: trunk/Ohana/src/relphot/src/select_images.c =================================================================== --- trunk/Ohana/src/relphot/src/select_images.c (revision 40019) +++ trunk/Ohana/src/relphot/src/select_images.c (revision 40291) @@ -251,20 +251,4 @@ image[nimage] = timage[i]; inSubset[i] = TRUE; - /* always allow 'few' images to succeed, if possible (new images / detections may have - * been added) */ - if (image[nimage].flags & ID_IMAGE_PHOTOM_FEW) { - image[nimage].flags &= ~(ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR); - } - if (RESET) { - if (RESET_ZEROPTS) { - if (!KEEP_UBERCAL || !(image[nimage].flags & ID_IMAGE_PHOTOM_UBERCAL)) { - image[nimage].Mcal = 0.0; - image[nimage].dMcal = NAN; - image[nimage].flags &= ~ID_IMAGE_PHOTOM_UBERCAL; - } - } - image[nimage].flags &= ~ID_IMAGE_PHOTOM_POOR; - image[nimage].ubercalDist = 1000; - } line_number[nimage] = i; nimage ++; Index: trunk/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 40019) +++ trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 40291) @@ -2,6 +2,6 @@ # if (0) -# define TEST_OBJ_ID 0x0000000e -# define TEST_CAT_ID 0x000076ee +# define TEST_OBJ_ID 0x000098f9 +# define TEST_CAT_ID 0x00005ad3 # else # define TEST_OBJ_ID 0 @@ -102,4 +102,8 @@ // 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 + +// Although I calculate McalAPER for exposures, I am only using McalPSF for chips. Note +// in StarOps.c:setMcalOutput I am setting measure->McalAPER to image->McalPSF for chips +// and warps, but not stacks // set mean of chip measurements (selected by photcode range for now): @@ -224,9 +228,8 @@ // overlaps). Msys + measure.Mcal is our best guess of the true magnitude Mmos = Mgrid = 0; - Mcal = measureT[k].Mcal; // check that this is zero for loaded REF value + Mcal = measureT[k].McalPSF; // 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); @@ -268,18 +271,21 @@ } - float Mkron, dMkron; + // NOTE: we need to calculate the averge chip Kron on each pass to be able to calibrate the stacks + float Mkron = PhotCatTiny (&measureT[k], MAG_CLASS_KRON); + float dMkron; if (isSetMrelFinal) { - Mkron = PhotCat (&measure[k], MAG_CLASS_KRON); dMkron = MAX (hypot(measure[k].dMkron, code->photomErrSys), MIN_ERROR); - if (CHECK_VALID_MAG(Mkron, dMkron)) { - int Nkron = results->kronData[Nsec].Nlist; - results->kronData[Nsec].flxlist[Nkron] = Mkron - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel - results->kronData[Nsec].errlist[Nkron] = dMkron; - results->kronData[Nsec].wgtlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0; - results->kronData[Nsec].ranking[Nkron] = measureRank[k]; - results->kronData[Nsec].measSeq[Nkron] = k; - results->kronData[Nsec].msklist[Nkron] = 0; - results->kronData[Nsec].Nlist ++; - } + } else { + dMkron = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR); + } + if (CHECK_VALID_MAG(Mkron, dMkron)) { + int Nkron = results->kronData[Nsec].Nlist; + results->kronData[Nsec].flxlist[Nkron] = Mkron - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel + results->kronData[Nsec].errlist[Nkron] = dMkron; + results->kronData[Nsec].wgtlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0; + results->kronData[Nsec].ranking[Nkron] = measureRank ? measureRank[k] : 0; + results->kronData[Nsec].measSeq[Nkron] = k; + results->kronData[Nsec].msklist[Nkron] = 0; + results->kronData[Nsec].Nlist ++; } @@ -377,8 +383,8 @@ } } - secfilt[Nsec].M = Mpsf + ZP; - secfilt[Nsec].dM = 0.6; - secfilt[Nsec].Mchisq = 0.0; - secfilt[Nsec].flags |= ID_SECF_USE_SYNTH; + secfilt[Nsec].MpsfChp = Mpsf + ZP; + secfilt[Nsec].dMpsfChp = 0.6; + secfilt[Nsec].Mchisq = 0.0; + secfilt[Nsec].flags |= ID_SECF_USE_SYNTH; continue; } @@ -386,13 +392,12 @@ // if too few valid measurements meet the minimum criteria, go to the next entry StatType *psfstats = &results->psfstats; - int Nranking = MAG_STATS_BY_RANKING (&results->psfData[Nsec], psfstats); - if (Nranking < Nminmeas) { + int NrankingPSF = MAG_STATS_BY_RANKING (&results->psfData[Nsec], psfstats); + if (NrankingPSF < Nminmeas) { secfilt[Nsec].flags |= ID_OBJ_FEW; } else { - secfilt[Nsec].M = psfstats->mean; - secfilt[Nsec].dM = psfstats->error; - secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN; - } - int minRankPSF = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10; + secfilt[Nsec].MpsfChp = psfstats->mean; + secfilt[Nsec].dMpsfChp = psfstats->error; + secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN; + } // when running -averages, we have no information about the images, so we cannot set this @@ -401,4 +406,13 @@ } + StatType *kronstats = &results->kronstats; + int NrankingKRON = MAG_STATS_BY_RANKING (&results->kronData[Nsec], kronstats); + if (NrankingKRON) { + secfilt[Nsec].MkronChp = kronstats->mean; + secfilt[Nsec].dMkronChp = kronstats->error; + secfilt[Nsec].sMkronChp = kronstats->sigma; + secfilt[Nsec].NusedKron = NrankingKRON; + } + if (isSetMrelFinal) { if ((average[0].objID == TEST_OBJ_ID) && (average[0].catID == TEST_CAT_ID)) { @@ -406,12 +420,9 @@ } - // mark the measurements matching this ranking - markMeasureByRanking (&results->psfData[Nsec], measure, minRankPSF, ID_MEAS_PHOTOM_PSF, ID_MEAS_MASKED_PSF); - - if (Nranking) { - secfilt[Nsec].Mstdev = psfstats->sigma; - secfilt[Nsec].Nused = psfstats->Nmeas; - secfilt[Nsec].Mmax = psfstats->max; - secfilt[Nsec].Mmin = psfstats->min; + if (NrankingPSF) { + secfilt[Nsec].sMpsfChp = psfstats->sigma; + secfilt[Nsec].Nused = psfstats->Nmeas; + secfilt[Nsec].Mmax = psfstats->max; + secfilt[Nsec].Mmin = psfstats->min; } secfilt[Nsec].psfQfMax = results->psfQfMax[Nsec]; @@ -421,25 +432,21 @@ MAX_NOT_NAN (psfQfPerfMax, secfilt[Nsec].psfQfPerfMax); + // mark the measurements matching this ranking + int minRankPSF = (NrankingPSF > 0) ? results->psfData[Nsec].ranking[0] : 10; + markMeasureByRanking (&results->psfData[Nsec], measure, minRankPSF, ID_MEAS_PHOTOM_PSF, ID_MEAS_MASKED_PSF); + + int minRankKron = (NrankingKRON > 0) ? results->kronData[Nsec].ranking[0] : 10; + markMeasureByRanking (&results->kronData[Nsec], measure, minRankKron, ID_MEAS_PHOTOM_KRON, ID_MEAS_MASKED_KRON); + StatType *apstats = &results->apstats; - Nranking = MAG_STATS_BY_RANKING (&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; - } - int minRankAper = (Nranking > 0) ? results->aperData[Nsec].ranking[0] : 10; + int NrankingAPER = MAG_STATS_BY_RANKING (&results->aperData[Nsec], apstats); + if (NrankingAPER) { + secfilt[Nsec].MapChp = apstats->mean; + secfilt[Nsec].dMapChp = apstats->error; + secfilt[Nsec].sMapChp = apstats->sigma; + secfilt[Nsec].NusedAp = NrankingAPER; + } + int minRankAper = (NrankingAPER > 0) ? results->aperData[Nsec].ranking[0] : 10; markMeasureByRanking (&results->aperData[Nsec], measure, minRankAper, ID_MEAS_PHOTOM_APER, ID_MEAS_MASKED_APER); - - StatType *kronstats = &results->kronstats; - Nranking = MAG_STATS_BY_RANKING (&results->kronData[Nsec], kronstats); - if (Nranking) { - secfilt[Nsec].Mkron = kronstats->mean; - secfilt[Nsec].dMkron = kronstats->error; - secfilt[Nsec].sMkron = kronstats->sigma; - secfilt[Nsec].NusedKron = Nranking; - } - int minRankKron = (Nranking > 0) ? results->kronData[Nsec].ranking[0] : 10; - markMeasureByRanking (&results->kronData[Nsec], measure, minRankKron, ID_MEAS_PHOTOM_KRON, ID_MEAS_MASKED_KRON); // does this object appear extended in > 50% of measurements? @@ -518,4 +525,9 @@ // 2) select the BEST detections per filter (regardless of PRIMARY) // 3) apply the zero point and AB->Jy transformations + +// I calculate McalAPER and McalPSF independently for stacks. I use McalAPER for Mkron +// and Map, and McalPSF for Mpsf. Note in StarOps.c:setMcalOutput I am setting +// measure->McalAPER to image->McalPSF for chips and warps, but not stacks + int setMrelAverageStack (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr) { @@ -528,7 +540,5 @@ SecFilt *secfilt = &catalog[0].secfilt[ave*Nsecfilt]; - off_t k; - - float Mcal = 0, Mmos = 0, Mgrid = 0, Finst = 0; + float McalPSF = 0, McalAPER = 0, Mmos = 0, Mgrid = 0; // set the primary projection cell and skycell for this coordinate @@ -546,6 +556,5 @@ int haveStackObject = FALSE; - int Ns; - for (Ns = 0; Ns < Nphotcodes; Ns++) { + for (int Ns = 0; Ns < Nphotcodes; Ns++) { int thisCode = photcodes[Ns][0].code; @@ -556,20 +565,37 @@ int haveStack = FALSE; - float psfQFbest = 0.0; - - off_t stackBestMeasure = -1; - off_t stackPrimaryMeasure = -1; - int isBad = FALSE; int isSuspect = FALSE; - int Nstack = 0; // number for this photcode - int NstackDet = 0; // number for this photcode - + int NstackMeas = 0; // number of stack measurements for this photcode + int NstackDet = 0; // number of stack detections for this photcode (not forced) + + int primaryLevelMax = 0; + off_t primaryEntryMax = -1; + float primaryValueMax = 0.0; + + int bestLevelMax = 0; + off_t bestEntryMax = -1; + float bestValueMax = 0.0; + + int Nprimary = 0; + + // reset all stack-related values for this secfilt: secfilt[Nsec].stackBestOff = -1; secfilt[Nsec].stackPrmryOff = -1; - - off_t meas = measureOffset; - for (k = 0; k < Nmeasure; k++, meas++) { + secfilt[Nsec].FpsfStk = NAN; + secfilt[Nsec].dFpsfStk = NAN; + secfilt[Nsec].FkronStk = NAN; + secfilt[Nsec].dFkronStk = NAN; + secfilt[Nsec].FapStk = NAN; + secfilt[Nsec].dFapStk = NAN; + secfilt[Nsec].MpsfStk = NAN; + secfilt[Nsec].MkronStk = NAN; + secfilt[Nsec].MapStk = NAN; + secfilt[Nsec].Nstack = 0; + secfilt[Nsec].NstackDet = 0; + secfilt[Nsec].flags &= ~ID_SECF_STACK_FLAGS; // reset the stack flags + + for (off_t k = 0; k < Nmeasure; k++) { // only examine gpc1 stack data @@ -581,6 +607,6 @@ if (code->equiv != thisCode) { continue; } - Nstack ++; - if ((measure[k].photFlags2 & 0x00000004) == 0) NstackDet ++; + NstackMeas ++; + if ((measure[k].photFlags2 & 0x00000004) == 0) NstackDet ++; // 0x04 = PM_SOURCE_MODE2_MATCHED // clear this bit for all measurements @@ -590,87 +616,101 @@ haveStack = TRUE; haveStackObject = TRUE; + + // if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad); + + int isPrimary = FALSE; - // ** find the PRIMARY measurement - - // if we request the primary (USE_TREE_FOR_PRIMARY), this is true if the measurement is from the - // primary skycell for this position - if (MatchImageSkycellID (meas, cat, tessID, projectID, skycellID)) { - stackPrimaryMeasure = k; + // ** is this a PRIMARY measurement? (there may be more than one) + + // if we request the primary (USE_TREE_FOR_PRIMARY), this is true if the measurement + // is from the primary skycell for this position. (note that MatchImageSkycellID + // requires the entry in the full Measure table, not just this object) + if (MatchImageSkycellID (measureOffset + k, cat, tessID, projectID, skycellID)) { measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY; - secfilt[Nsec].stackPrmryOff = meas; - myAssert (secfilt[Nsec].stackPrmryOff <= catalog[0].Nmeasure, "stackPrmryOff out of range"); - } - - // ** now choose the BEST measurements (may also be PRIMARY) - - // ensure that we at least have a single best measure - if (stackBestMeasure == -1) stackBestMeasure = k; - - // choose the best psfQFperf value for the BEST measurement - if (isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > psfQFbest)) { - psfQFbest = measure[k].psfQFperf; - stackBestMeasure = k; - } - // ... UNLESS psfQFperf > 0.98 for the primary, in which case just use the primary. - if ((measure[k].dbFlags & ID_MEAS_STACK_PRIMARY) && isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > 0.98)) { - psfQFbest = 1000; // force this to be the best entry - stackBestMeasure = k; - } - - if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad); - + isPrimary = TRUE; + Nprimary ++; + } + + // soften the error floor (can dM be 0.0?) + int SNvalue = isfinite(measure[k].dM) ? 1.0 / hypot (measure[k].dM, MIN_ERROR) : NAN; + int psfQFperfAboveLimit = isfinite(SNvalue) && isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > 0.95); + + // ** determine the BEST level + int bestLevel = 0; + if (!isPrimary && !psfQFperfAboveLimit) bestLevel = 1; + if ( isPrimary && !psfQFperfAboveLimit) bestLevel = 2; + if (!isPrimary && psfQFperfAboveLimit) bestLevel = 3; + if ( isPrimary && psfQFperfAboveLimit) bestLevel = 4; + + // here is the rule for choosing the best value: + float bestValue = (bestLevel > 2) ? SNvalue : measure[k].psfQFperf; + + // if we have not reached this level before, set the new level + if (bestLevelMax < bestLevel) { + bestValueMax = bestValue; + bestLevelMax = bestLevel; + bestEntryMax = k; + } + // if we have reached this level before, set the new value if we beat the old one + if ((bestLevelMax == bestLevel) && (bestValueMax < bestValue)) { + bestValueMax = bestValue; + bestEntryMax = k; + } + myAssert (bestEntryMax > -1, "this should not happen"); + + // ** determine the PRIMARY level + int primaryLevel = 0; + if ( isPrimary && !psfQFperfAboveLimit) primaryLevel = 1; + if ( isPrimary && psfQFperfAboveLimit) primaryLevel = 2; + + // here is the rule for choosing the PRIMARY value: + float primaryValue = 0.0; + if (primaryLevel == 1) { primaryValue = measure[k].psfQFperf; } + if (primaryLevel == 2) { primaryValue = SNvalue; } + + // if we have not reached this level before, set the new level + if (isPrimary && (primaryLevelMax < primaryLevel)) { + primaryValueMax = primaryValue; + primaryLevelMax = primaryLevel; + primaryEntryMax = k; + } + // if we have reached this level before, set the new value if we beat the old one + if (isPrimary && (primaryLevelMax == primaryLevel) && (primaryValueMax < primaryValue)) { + primaryValueMax = primaryValue; + primaryEntryMax = k; + } + } + + // if we do not have a stack measurement for this photcode, skip everything below + if (!haveStack) continue; + + // now we have a BEST and a PRIMARY entry (may be the same entry) + + off_t meas = bestEntryMax; + off_t measSeq = meas + measureOffset; + + myAssert (meas >= 0, "this should not happen"); + + // measurements without an image are either external reference photometry or + // data for which the associated image has not been loaded (probably because of + // overlaps). we only want measurements associated with stack images in this loop + + // match measurement to its image (this is just a check, right?) + if (getImageEntry (measSeq, 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 - // overlaps). we only want measurements associated with stack images in this loop - - // match measurement to its image (this is just a check, right?) - 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 - // 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 - } else { - Mcal = getMcal_alt (meas, cat, flatcorr, measure[k].Xccd, measure[k].Yccd); - if (isnan(Mcal)) SKIP_THIS_MEAS_STACK(Ncal); - Mmos = getMmos (meas, cat); - if (isnan(Mmos)) SKIP_THIS_MEAS_STACK(Nmos); - Mgrid = getMgrid (meas, cat); - if (isnan(Mgrid)) SKIP_THIS_MEAS_STACK(Ngrid); - } - - // NOTE: negative and insignificant flux values are allowed, but not NAN flux values - Finst = PhotFluxInst (&measure[k], MAG_CLASS_PSF); - if (isnan(Finst)) SKIP_THIS_MEAS_STACK(Ninst); - - // data quality assessment - isBad |= (measure[k].photFlags & code->photomBadMask); - isBad |= (measure[k].psfQF < 0.85); - isBad |= isnan(measure[k].psfQF); - isBad |= measure[k].dM > 0.2; // S/N < 5.0 - - isSuspect |= (measure[k].photFlags & code->photomPoorMask); - isSuspect |= (measure[k].psfQFperf < 0.85); - } - - if (!haveStack) continue; - - if (!isSuspect && !isBad) { - NstackGood ++; - } - if (isSuspect && !isBad) { - NstackSuspect ++; - } - - // measurements which are bad will not have a valid stack entry and are skipped - k = (stackBestMeasure >= 0) ? stackBestMeasure : stackPrimaryMeasure; - if (k < 0) continue; - - // we are now populating stackDetectID not stack Image ID in secfilt - // ID = (stackPrimaryMeasureMin >= 0) ? stackPrimaryIDmin : stackCenterIDmin; - // ID = measure[k].extID; // for the stack, this is the stackDetectID - - // get the zero point for the selected image - float zp = PhotZeroPoint (&measure[k], &average[0], &secfilt[0]) - (Mcal + Mmos + Mgrid); + // overlaps). Msys + measure.Mcal is our best guess of the true magnitude + McalPSF = measure[meas].McalPSF; // check that this is zero for loaded REF value + McalAPER = McalPSF; // check that this is zero for loaded REF value + Mmos = 0.0; + Mgrid = 0.0; + } else { + McalPSF = getMcal (measSeq, cat, MAG_CLASS_PSF); + McalAPER = USE_MCAL_PSF_FOR_STACK_APER ? getMcal (measSeq, cat, MAG_CLASS_PSF) : getMcal (measSeq, cat, MAG_CLASS_KRON); + Mmos = getMmos (measSeq, cat); + Mgrid = getMgrid (measSeq, cat); + } + + // ** Here is the math to relate mag,zp to flux // flux_cgs : erg sec^1 cm^-2 Hz^-1 @@ -703,46 +743,81 @@ // flux_Jy = flux_inst * ten(-0.4*ZP + 3.56) + // get the zero point for the selected image + // Use a different zero point for the PSF-like and APERTURE-like magnitudes + float zpPSF = PhotZeroPoint (&measure[meas], &average[0], &secfilt[0]) - (McalPSF + Mmos + Mgrid); + float zpAPER = PhotZeroPoint (&measure[meas], &average[0], &secfilt[0]) - (McalAPER + Mmos + Mgrid); + // 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[meas].FluxPSF; + secfilt[Nsec].dFpsfStk = zpFactorPSF * measure[meas].dFluxPSF; + secfilt[Nsec].FkronStk = zpFactorAPER * measure[meas].FluxKron; + secfilt[Nsec].dFkronStk = zpFactorAPER * measure[meas].dFluxKron; + secfilt[Nsec].FapStk = zpFactorAPER * measure[meas].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[meas].dFluxAp; + secfilt[Nsec].dFapStk = zpFactorAPER * measure[meas].dFluxPSF; // Jy to AB mags - secfilt[Nsec].MpsfStk = (measure[k].FluxPSF > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FpsfStk) : NAN; - secfilt[Nsec].MkronStk = (measure[k].FluxKron > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FkronStk) : NAN; - secfilt[Nsec].MapStk = (measure[k].FluxAp > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FapStk) : NAN; - - secfilt[Nsec].stackBestOff = k + measureOffset; + secfilt[Nsec].MpsfStk = (measure[meas].FluxPSF > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FpsfStk) : NAN; + secfilt[Nsec].MkronStk = (measure[meas].FluxKron > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FkronStk) : NAN; + secfilt[Nsec].MapStk = (measure[meas].FluxAp > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FapStk) : NAN; + + // record the measurement which gave the best value + secfilt[Nsec].stackBestOff = (bestEntryMax == -1) ? -1 : bestEntryMax + measureOffset; myAssert (secfilt[Nsec].stackBestOff <= catalog[0].Nmeasure, "stackBestOff out of range"); - secfilt[Nsec].Nstack = Nstack; + // record the selected primary measurement + secfilt[Nsec].stackPrmryOff = (primaryEntryMax == -1) ? -1 : primaryEntryMax + measureOffset; + myAssert (secfilt[Nsec].stackPrmryOff <= catalog[0].Nmeasure, "stackPrmryOff out of range"); + + // this is the selected measurement used by secfilt[] + measure[meas].dbFlags |= ID_MEAS_STACK_PHOT_SRC; + if (Nprimary) { + secfilt[Nsec].flags |= ID_SECF_STACK_PRIMARY; + if (Nprimary > 1) { + secfilt[Nsec].flags |= ID_SECF_STACK_PRIMARY_MULTIPLE; + } + } + + // stack measurement used for 'best' was a detection (not forced from the other bands) + if ((measure[bestEntryMax].photFlags2 & 0x00000004) == 0) { + secfilt[Nsec].flags |= ID_SECF_STACK_BESTDET; + } + + // stack measurement used for 'primary' was a detection (not forced from the other bands) + if ((primaryEntryMax >= 0) && ((measure[primaryEntryMax].photFlags2 & 0x00000004) == 0)) { + secfilt[Nsec].flags |= ID_SECF_STACK_PRIMDET; + } + + secfilt[Nsec].flags |= ID_SECF_HAS_PS1_STACK; + + // NOTE: negative and insignificant flux values are allowed, but not NAN flux values + // float Finst = PhotFluxInst (&measure[meas], MAG_CLASS_PSF); + + // data quality assessment + isBad |= (measure[meas].photFlags & photcodes[Ns][0].photomBadMask); + isBad |= (measure[meas].psfQF < 0.85); + isBad |= isnan(measure[meas].psfQF); + isBad |= measure[meas].dM > 0.2; // S/N < 5.0 + + isSuspect |= (measure[meas].photFlags & photcodes[Ns][0].photomPoorMask); + isSuspect |= (measure[meas].psfQFperf < 0.85); + + if (!isSuspect && !isBad) { + NstackGood ++; + } + if (isSuspect && !isBad) { + NstackSuspect ++; + } + + secfilt[Nsec].Nstack = NstackMeas; secfilt[Nsec].NstackDet = NstackDet; - // this is the measurement used by secfilt[] - measure[k].dbFlags |= ID_MEAS_STACK_PHOT_SRC; - if (k == stackPrimaryMeasure) { - secfilt[Nsec].flags |= ID_SECF_STACK_PRIMARY; - } - - // stack measurement used for 'best' was a detection (not forced from the other bands) - if ((measure[k].photFlags2 & 0x00000004) == 0) { - secfilt[Nsec].flags |= ID_SECF_STACK_BESTDET; - } - // stack measurement used for 'primary' was a detection (not forced from the other bands) - if ((stackPrimaryMeasure >= 0) && ((measure[stackPrimaryMeasure].photFlags2 & 0x00000004) == 0)) { - secfilt[Nsec].flags |= ID_SECF_STACK_PRIMDET; - } - - secfilt[Nsec].flags |= ID_SECF_HAS_PS1_STACK; } // Nsecfilt loop @@ -755,5 +830,5 @@ int PrimaryIsBest = TRUE; - for (Ns = 0; Ns < Nphotcodes; Ns++) { + for (int Ns = 0; Ns < Nphotcodes; Ns++) { if (!(secfilt[Ns].flags & ID_SECF_HAS_PS1_STACK)) continue; // no stack detection in PS1, nothing is best if (secfilt[Ns].flags & ID_SECF_STACK_PRIMARY) continue; // primary stack detection is best @@ -780,4 +855,9 @@ // * no grid, no mosaic, no 2MASS, no SYNTH, no Ubercal, no flatcorr // analysis is done on flux, not mags (as the faintest objects will be nearly insignificant) + +// Although I calculate McalAPER for exposures, I am only using McalPSF for warps.. +// Note in StarOps.c:setMcalOutput I am setting measure->McalAPER to image->McalPSF for +// chips and warps, but not stacks + int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results) { OHANA_UNUSED_PARAM(flatcorr); @@ -857,9 +937,7 @@ // data for which the associated image has not been loaded (probably because of // overlaps). Msys + measure.Mcal is our best guess of the true magnitude - Mcal = measure[k].Mcal; // check that this is zero for loaded REF value + Mcal = measure[k].McalPSF; // 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: trunk/Ohana/src/relphot/src/setMrelFinal.c =================================================================== --- trunk/Ohana/src/relphot/src/setMrelFinal.c (revision 40019) +++ trunk/Ohana/src/relphot/src/setMrelFinal.c (revision 40291) @@ -65,10 +65,12 @@ if (catalog[0].measure[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) { if (!KEEP_UBERCAL) { - catalog[0].measure[m].Mcal = 0.0; + catalog[0].measure[m].McalPSF = 0.0; + catalog[0].measure[m].McalAPER = 0.0; catalog[0].measure[m].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL; } } else { if (RESET_ZEROPTS && (getImageEntry (m, 0) >= 0)) { - catalog[0].measure[m].Mcal = 0.0; + catalog[0].measure[m].McalPSF = 0.0; + catalog[0].measure[m].McalAPER = 0.0; } } @@ -105,5 +107,5 @@ ALLOCATE (catalog[0].measureRank, char, catalog[0].Nmeasure); setMeasureRank (catalog); - setMrelOutput (catalog, 1, flatcorr); // sets the values secfilt.M = + setMrelOutput (catalog, 1, flatcorr); // sets the values secfilt.MpsfChp = setMcalOutput (catalog, 1, flatcorr); // sets measure.Mcal = image.Mcal @@ -226,5 +228,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; } Index: trunk/Ohana/src/relphot/src/share_image_mags.c =================================================================== --- trunk/Ohana/src/relphot/src/share_image_mags.c (revision 40019) +++ trunk/Ohana/src/relphot/src/share_image_mags.c (revision 40291) @@ -74,8 +74,9 @@ continue; } - images[seq].Mcal = image_mags[i].Mcal; + images[seq].McalPSF = image_mags[i].McalPSF; + images[seq].McalAPER = image_mags[i].McalAPER; images[seq].dMcal = image_mags[i].dMcal; images[seq].dMagSys = image_mags[i].dMagSys; - images[seq].Xm = image_mags[i].Xm; + images[seq].McalChiSq = image_mags[i].McalChiSq; images[seq].nFitPhotom = image_mags[i].nFitPhotom; images[seq].flags = image_mags[i].flags; @@ -91,8 +92,9 @@ int set_image_mags (ImageMag *image_mags, Image *image) { - image_mags->Mcal = image->Mcal; + image_mags->McalPSF = image->McalPSF; + image_mags->McalAPER = image->McalAPER; image_mags->dMcal = image->dMcal; image_mags->dMagSys = image->dMagSys; - image_mags->Xm = image->Xm; + image_mags->McalChiSq = image->McalChiSq; image_mags->nFitPhotom = image->nFitPhotom; image_mags->flags = image->flags; Index: trunk/Ohana/src/relphot/src/share_mean_mags.c =================================================================== --- trunk/Ohana/src/relphot/src/share_mean_mags.c (revision 40019) +++ trunk/Ohana/src/relphot/src/share_mean_mags.c (revision 40291) @@ -111,5 +111,5 @@ if (Nsec < 0) continue; - catalog[catSeq].secfilt[objSeq*Nsecfilt + Nsec].M = meanmags[i].M; + catalog[catSeq].secfilt[objSeq*Nsecfilt + Nsec].MpsfChp = meanmags[i].M; } free (meanmags); @@ -122,6 +122,6 @@ int set_mean_mags (MeanMag *meanmags, AverageTiny *average, SecFilt *secfilt, int Nsec) { - meanmags->M = secfilt->M; - meanmags->dM = secfilt->dM; + meanmags->M = secfilt->MpsfChp; + meanmags->dM = secfilt->dMpsfChp; meanmags->Mchisq = secfilt->Mchisq; meanmags->Nsec = Nsec; // key to secfilt entry Index: trunk/Ohana/src/relphot/src/synthetic_mags.c =================================================================== --- trunk/Ohana/src/relphot/src/synthetic_mags.c (revision 40019) +++ trunk/Ohana/src/relphot/src/synthetic_mags.c (revision 40291) @@ -37,6 +37,6 @@ // (r - i > 0.5) : w = r + 0.268 - 0.435 (r-i) - 0.078(r-i)^2 - float Mr = secfilt[NSr].M; - float Mi = secfilt[NSi].M; + float Mr = secfilt[NSr].MpsfChp; + float Mi = secfilt[NSi].MpsfChp; if (!isfinite(Mr)) return FALSE;