Index: /trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- /trunk/Ohana/src/relphot/include/relphot.h (revision 41556) +++ /trunk/Ohana/src/relphot/include/relphot.h (revision 41557) @@ -598,4 +598,5 @@ void FitDataSetAlloc PROTO((FitDataSet *dataset, int Nmax, int order, int Nbootstrap)); void FitDataSetAddPriors PROTO((FitDataSet *dataset)); +StatType FitDataSetSoften PROTO((FitDataSet *dataset, int Nvalues)); unsigned int *ReadTGroupFile PROTO((FILE *f, int *nelem)); Index: /trunk/Ohana/src/relphot/src/ImageOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/ImageOps.c (revision 41556) +++ /trunk/Ohana/src/relphot/src/ImageOps.c (revision 41557) @@ -784,8 +784,6 @@ } - // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not) - // XXX: apply airmass from above and fit only zp? - - // XXX soften the stats? + // soften the errors based on the scatter + FitDataSetSoften (&psfStars, Nref); fit1d_irls (&psfStars, Nref); Index: /trunk/Ohana/src/relphot/src/MagResidSave.c =================================================================== --- /trunk/Ohana/src/relphot/src/MagResidSave.c (revision 41556) +++ /trunk/Ohana/src/relphot/src/MagResidSave.c (revision 41557) @@ -12,4 +12,6 @@ float dMresid; float Mchisq; + float Mmedian; + float Maltsig; int Nfit; } MagResidResult; @@ -39,8 +41,9 @@ } - // we are making a 0-order fit and not doing bootstrap analysis, using OLS + // we are making a 0-order fit and not doing bootstrap analysis FitDataSet psfStars; FitDataSetAlloc (&psfStars, Nmax, 0, 0); - psfStars.MaxIterations = 0; + // By default, we use IRLS; if MaxIterations = 0, then we use OLS + // psfStars.MaxIterations = 0; gfits_init_header (&header); @@ -56,4 +59,6 @@ gfits_define_bintable_column (&theader, "E", "MAG_ERROR", "magnitude error", "magnitudes", 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "MAG_CHISQ", "magnitude resid chisq", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_MEDIAN", "median magnitude residual", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_ALTSIG", "68 percentile range as sigma", "magnitudes", 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "AIRMASS", "", "unitless", 1.0, 0.0); gfits_define_bintable_column (&theader, "D", "MJD", "", "unitless", 1.0, 0.0); @@ -62,4 +67,5 @@ gfits_define_bintable_column (&theader, "J", "FLAGS", "image flags", "unitless", 1.0, 0.0); gfits_define_bintable_column (&theader, "J", "NFIT", "number of stars fitted", NULL, 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "NMEAS", "number of stars on image", NULL, 1.0, FT_BZERO_INT32); // generate the output array that carries the data @@ -73,4 +79,6 @@ ALLOCATE_PTR (dMresid, float, Nimage); ALLOCATE_PTR (Mchisq, float, Nimage); + ALLOCATE_PTR (Mmedian, float, Nimage); + ALLOCATE_PTR (Maltsig, float, Nimage); ALLOCATE_PTR (airmass, float, Nimage); ALLOCATE_PTR (mjd, double, Nimage); @@ -79,4 +87,5 @@ ALLOCATE_PTR (flags, int, Nimage); ALLOCATE_PTR (Nfit, int, Nimage); + ALLOCATE_PTR (Nmeas, int, Nimage); // XXXX fill in the vectors here @@ -94,4 +103,6 @@ Mchisq[i] = NAN; dMresid[i] = NAN; + Mmedian[i] = NAN; + Maltsig[i] = NAN; Nfit[i] = 0; Nmeas[i] = N_onImage[i]; @@ -110,13 +121,17 @@ Mchisq[i] = result.Mchisq; dMresid[i] = result.dMresid; + Mmedian[i] = result.Mmedian; + Maltsig[i] = result.Maltsig; Nfit[i] = result.Nfit; } // add the columns to the output array - gfits_set_bintable_column (&theader, &ftable, "RA", Rs, Nimage); - gfits_set_bintable_column (&theader, &ftable, "DEC", Ds, Nimage); + gfits_set_bintable_column (&theader, &ftable, "RA", Rs, Nimage); + gfits_set_bintable_column (&theader, &ftable, "DEC", Ds, Nimage); gfits_set_bintable_column (&theader, &ftable, "MAG_RESID", Mresid, Nimage); gfits_set_bintable_column (&theader, &ftable, "MAG_ERROR", dMresid, Nimage); gfits_set_bintable_column (&theader, &ftable, "MAG_CHISQ", Mchisq, Nimage); + gfits_set_bintable_column (&theader, &ftable, "MAG_MEDIAN", Mmedian, Nimage); + gfits_set_bintable_column (&theader, &ftable, "MAG_ALTSIG", Maltsig, Nimage); gfits_set_bintable_column (&theader, &ftable, "AIRMASS", airmass, Nimage); gfits_set_bintable_column (&theader, &ftable, "MJD", mjd, Nimage); @@ -125,4 +140,5 @@ gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags, Nimage); gfits_set_bintable_column (&theader, &ftable, "NFIT", Nfit, Nimage); + gfits_set_bintable_column (&theader, &ftable, "NMEAS", Nmeas, Nimage); free (Rs); @@ -131,4 +147,6 @@ free (dMresid ); free (Mchisq ); + free (Mmedian ); + free (Maltsig ); free (airmass ); free (mjd ); @@ -137,4 +155,5 @@ free (flags ); free (Nfit ); + free (Nmeas ); FILE *f = fopen (filename, "w"); @@ -185,4 +204,6 @@ result.dMresid = NAN; result.Mchisq = NAN; + result.Mmedian = NAN; + result.Maltsig = NAN; result.Nfit = 0; @@ -239,7 +260,12 @@ result.dMresid = NAN; result.Mchisq = NAN; + result.Mmedian = NAN; + result.Maltsig = NAN; result.Nfit = Nref; return result; } + + StatType stats = FitDataSetSoften (psfStars, Nref); + double altSigma = (stats.Upper80 - stats.Lower20) / 1.6; // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not) @@ -249,4 +275,6 @@ result.Nfit = psfStars->Nmeas; result.Mchisq = psfStars->chisq; + result.Mmedian = stats.median; + result.Maltsig = altSigma; return result; Index: /trunk/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/MosaicOps.c (revision 41556) +++ /trunk/Ohana/src/relphot/src/MosaicOps.c (revision 41557) @@ -1161,4 +1161,5 @@ // in GOOD_MOSAIC mode, we fit good mosaics ignoring night state + // This case should only be used when not fitting nights / tgroups if ((MOSAIC_ZPT_MODE == MOSAIC_ZPT_MODE_GOOD_MOSAIC)) { if (badMosaic) useMmos = FALSE; // do not fit bad mosaics @@ -1312,15 +1313,7 @@ } - StatType stats; - liststats_setmode (&stats, "MEDIAN"); - liststats (psfStars->alldata->yVector, NULL, NULL, N, &stats); - double altSigma = (stats.Upper80 - stats.Lower20) / 1.6; // 20% to 80% encompasses 60% of the values, corresponds to the range (-0.85 sigma : +0.85 sigma) - - // soften the individual errors with 10% of the scatter above - for (j = 0; j < N; j++) { - double newSigma = hypot(psfStars->alldata->dyVector[j], 0.1*altSigma); - psfStars->alldata->dyVector[j] = newSigma; - } - + // soften the errors based on the scatter + FitDataSetSoften (psfStars, N); + fit1d_irls (psfStars, N); fit1d_irls (brightStars, Nbright); Index: /trunk/Ohana/src/relphot/src/TGroupOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/TGroupOps.c (revision 41556) +++ /trunk/Ohana/src/relphot/src/TGroupOps.c (revision 41557) @@ -927,18 +927,9 @@ */ - // use liststats to find the 20-pct, median, 80-pct points - StatType stats; - liststats_setmode (&stats, "MEDIAN"); - liststats (psfStars->alldata->yVector, NULL, NULL, N, &stats); - double altSigma = (stats.Upper80 - stats.Lower20) / 1.6; // 20% to 80% encompasses 60% of the values, corresponds to the range (-0.85 sigma : +0.85 sigma) + StatType stats = FitDataSetSoften (psfStars, N); + double altSigma = (stats.Upper80 - stats.Lower20) / 1.6; TGTimes *mygroup = (TGTimes *) myTGroup->parent; fprintf (stderr, "TGroup Stats %f,%d : %d %d %6.3f %6.3f\n", ohana_sec_to_mjd(mygroup->start), myTGroup->photcode, (int) N, (int) myTGroup->Nmeasure, stats.median, altSigma); - - // soften the individual errors with 10% of the scatter above - for (j = 0; j < N; j++) { - double newSigma = hypot(psfStars->alldata->dyVector[j], 0.1*altSigma); - psfStars->alldata->dyVector[j] = newSigma; - } // do anything special with identified good nights? Index: /trunk/Ohana/src/relphot/src/bcatalog.c =================================================================== --- /trunk/Ohana/src/relphot/src/bcatalog.c (revision 41556) +++ /trunk/Ohana/src/relphot/src/bcatalog.c (revision 41557) @@ -106,6 +106,6 @@ // skip garbage measurements - if (isnan(catalog[0].measure[offset].psfQF) || (catalog[0].measure[offset].psfQF < 0.85)) { Npsfqf ++; continue; } - if (isnan(catalog[0].measure[offset].psfQFperf) || (catalog[0].measure[offset].psfQFperf < 0.85)) { Npsfqf ++; continue; } + if (isnan(catalog[0].measure[offset].psfQF) || (catalog[0].measure[offset].psfQF < 0.95)) { Npsfqf ++; continue; } + if (isnan(catalog[0].measure[offset].psfQFperf) || (catalog[0].measure[offset].psfQFperf < 0.95)) { Npsfqf ++; continue; } if (isnan(catalog[0].measure[offset].M)) { Nnan ++; continue; } Index: /trunk/Ohana/src/relphot/src/fit1d_irls.c =================================================================== --- /trunk/Ohana/src/relphot/src/fit1d_irls.c (revision 41556) +++ /trunk/Ohana/src/relphot/src/fit1d_irls.c (revision 41557) @@ -368,2 +368,20 @@ } } + +StatType FitDataSetSoften (FitDataSet *dataset, int Nvalues) { + + // use liststats to find the 20-pct, median, 80-pct points + StatType stats; + liststats_setmode (&stats, "MEDIAN"); + liststats (dataset->alldata->yVector, NULL, NULL, Nvalues, &stats); + + double altSigma = (stats.Upper80 - stats.Lower20) / 1.6; // 20% to 80% encompasses 60% of the values, corresponds to the range (-0.85 sigma : +0.85 sigma) + + // soften the individual errors with 10% of the scatter above + for (int j = 0; j < Nvalues; j++) { + double newSigma = hypot(dataset->alldata->dyVector[j], 0.1*altSigma); + dataset->alldata->dyVector[j] = newSigma; + } + + return stats; +} Index: /trunk/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- /trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 41556) +++ /trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 41557) @@ -267,4 +267,6 @@ if (USE_REFERENCE_WEIGHT && (code->type == PHOT_REF)) modifiedWeight = code->photomErrSys; // we are overloading this field for now + // XXX I should probably do something more clever with the measured scatter in the nights, mosaics, images + // combine the various errors in quadrature here: float dMsys = hypot (code->photomErrSys, dMgrp);