Index: /trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- /trunk/Ohana/src/relphot/include/relphot.h (revision 41554) +++ /trunk/Ohana/src/relphot/include/relphot.h (revision 41555) @@ -88,4 +88,7 @@ float dMcal; float dMsys; + float stdev; + float dMmin; + float dMmax; float McalChiSq; float secz; @@ -107,4 +110,6 @@ float dMsys; float stdev; + float dMmin; + float dMmax; float McalChiSq; unsigned int nFitPhotom; @@ -155,6 +160,8 @@ double min; double max; + double Upper90; double Upper80; double Lower20; + double Lower10; double total; int Nmeas; @@ -415,4 +422,5 @@ int TEST_IMAGE2; +int MANUAL_ITERATION; int NLOOP; int NGRID; @@ -643,4 +651,6 @@ int setMrelOutput PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); +void setMcalTest PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); + int setMave PROTO((Catalog *catalog, int Ncatalog)); void set_ZP PROTO((double ZERO)); @@ -835,4 +845,5 @@ void dump_tgroups (Catalog *catalog, int Npass); void dump_catalog (Catalog *catalog, off_t c, int Npass); +void dump_tgroup_imstats (int Npass); void SetZptIteration (int current); Index: /trunk/Ohana/src/relphot/src/ImageOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/ImageOps.c (revision 41554) +++ /trunk/Ohana/src/relphot/src/ImageOps.c (revision 41555) @@ -1,4 +1,3 @@ # include "relphot.h" -# define BASIC_STATS 0 // the MeasureToImage, ImageToCatalog, and ImageToMeasure arrays are a substantial part of the memory footprint. @@ -631,4 +630,7 @@ int isMosaic = isMosaicChip(image[i].photcode); + // unset at start by default (set if actually used below) + image[i].flags &= ~ID_IMAGE_IMAGE_PHOTCAL; + // if requested, freeze mosaic chips: if (FREEZE_IMAGES && isMosaic) continue; @@ -719,4 +721,6 @@ } + float Moff = Mmos + Mgrp + Mgrid - Mflat; + PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode); if (!code) goto skip; @@ -727,5 +731,5 @@ skip: - psfStars.alldata-> yVector[Nref] = MsysPSF - MrelPSF - Mmos - Mgrp - Mgrid + Mflat; + psfStars.alldata-> yVector[Nref] = MsysPSF - MrelPSF - Moff; psfStars.alldata->dyVector[Nref] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); @@ -734,5 +738,5 @@ if (isfinite(MrelKron) && isfinite(MsysKron)) { - kronStars.alldata-> yVector[Nkron] = MsysKron - MrelKron - Mmos - Mgrp - Mgrid + Mflat; + kronStars.alldata-> yVector[Nkron] = MsysKron - MrelKron - Moff; kronStars.alldata->dyVector[Nkron] = psfStars.alldata->dyVector[Nref]; Nkron ++; @@ -765,4 +769,5 @@ } } + if (mark) continue; // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not) @@ -776,6 +781,9 @@ // 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? - fit1d_irls (&kronStars, Nkron); - image[i].McalAPER = kronStars.bSaveArray[0][0]; + // XXX EAM : temporarily use this to track calibration quality + if (0) { + fit1d_irls (&kronStars, Nkron); + image[i].McalAPER = kronStars.bSaveArray[0][0]; + } if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) { @@ -795,4 +803,6 @@ image[i].McalPSF = 0.0; } + + image[i].flags |= ID_IMAGE_IMAGE_PHOTCAL; // set this flag // minUbercalDist calculated here is the min value for any star owned by this image @@ -814,4 +824,120 @@ FitDataSetFree (&psfStars); + return; +} + +/* determine McalTEST values for all images -- this is not used to set measure.Mcal, but only to test */ +void setMcalTest (Catalog *catalog, FlatCorrectionTable *flatcorr) { + + off_t i, j, m, c, n; + int Nfew, Nbad, Nmos, Ngrp, Nrel, Ngrid, Nsys; + + int Nsecfilt = GetPhotcodeNsecfilt (); + + off_t Nmax = 0; + for (i = 0; i < Nimage; i++) { + Nmax = MAX (Nmax, N_onImage[i]); + } + + // we are making a 0-order fit and not doing bootstrap analysis: + FitDataSet psfStars; + FitDataSetAlloc (&psfStars, Nmax, 0, 0); + + // for testing, need to allow this to be optional + if (FALSE && UseStandardOLS(ZPT_IMAGES)) { + psfStars.MaxIterations = 0; + } + psfStars.MaxIterations = 0; + + Nfew = Nbad = Nmos = Ngrp = Ngrid = Nrel = Nsys = 0; + + int Ncalibrated = 0; + for (i = 0; i < Nimage; i++) { + if (image[i].photcode == 0) continue; // skip the PHU images + + off_t Nref = 0; // number of stars used to measure McalPSF + + if (N_onImage[i] == 0) { + // fprintf (stderr, "image missing detections? %s %d "OFF_T_FMT"\n", image[i].name, image[i].nstar, N_onImage[i]); + continue; + } + + for (j = 0; j < N_onImage[i]; j++) { + + m = ImageToMeasure[i][j]; + c = ImageToCatalog[i][j]; + + if (catalog[c].measureT[m].dbFlags & MEAS_BAD) { + Nbad++; + continue; + } + float Mmos = getMmos (m, c); + if (isnan(Mmos)) { + Nmos ++; + continue; + } + float Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass, NULL); // ignore error for now? + if (isnan(Mgrp)) { + Ngrp ++; + continue; + } + float Mgrid = getMgrid (m, c); + if (isnan(Mgrid)) { + Ngrid++; + continue; + } + + // MrelPSF is the average magnitude for this star. + float MrelPSF = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); + if (isnan(MrelPSF)) { + Nrel ++; + continue; + } + + // image.Mcal is not supposed to include the flat-field correction, so we need to + // apply that offset as well here for this image (in other words, each detection is + // being compared to the model, excluding the zero point, Mcal. The model includes + // the flat-correction. NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat) + + float Mflat = getMflat (m, c, flatcorr, catalog); + + // get the PSF magnitude for thie measurement, with airmass slope applied + n = catalog[c].measureT[m].averef; + float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); + if (isnan(MsysPSF)) { + Nsys++; + continue; + } + + float Moff = Mmos + Mgrp + Mgrid - Mflat; + + psfStars.alldata-> yVector[Nref] = MsysPSF - MrelPSF - Moff; + psfStars.alldata->dyVector[Nref] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); + Nref++; + } + + if (Nref < 5) { + // fprintf (stderr, "very few detections? %s %d "OFF_T_FMT"\n", image[i].name, image[i].nstar, Nref); + // continue; + image[i].McalAPER = NAN; + image[i].dMcal = NAN; + image[i].nFitPhotom = Nref; + image[i].McalChiSq = NAN; + continue; + } + + // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not) + fit1d_irls (&psfStars, Nref); + image[i].McalAPER = psfStars.bSaveArray[0][0]; + image[i].dMcal = psfStars.bSigma[0]; + image[i].nFitPhotom = psfStars.Nmeas; + image[i].McalChiSq = psfStars.chisq; + Ncalibrated ++; + } + + fprintf (stderr, "%d images calibrated\n", Ncalibrated); + fprintf (stderr, "%d images marked having too few measurements (Nbad: %d, Nmos: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", Nfew, Nbad, Nmos, Ngrid, Nrel, Nsys); + + FitDataSetFree (&psfStars); return; } Index: /trunk/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/MosaicOps.c (revision 41554) +++ /trunk/Ohana/src/relphot/src/MosaicOps.c (revision 41555) @@ -685,4 +685,12 @@ // if it helps, note that ubercal uses a single zp per exposure, so the mean of those values is the same as the value + if (!isfinite(image[m].McalPSF)) { + image[m].McalPSF = 0.0; + image[m].McalAPER = 0.0; + image[m].dMcal = 0.0; + image[m].McalChiSq = 0.0; + fprintf (stderr, "warning: resetting NAN value for Mcal %s\n", image[m].name); + } + McalPSF += image[m].McalPSF; McalAPER += image[m].McalAPER; @@ -694,4 +702,5 @@ image[m].dMcal = NAN; image[m].McalChiSq = NAN; + } dS /= MosaicN_Image[i]; @@ -734,14 +743,19 @@ for (j = 0; j < MosaicN_Image[i]; j++) { im = MosaicToImage[i][j]; - 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; + if (mosaic[i].flags & ID_IMAGE_MOSAIC_PHOTCAL) { + 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); image[im].flags |= (mosaic[i].flags & ID_IMAGE_MOSAIC_POOR); + image[im].flags |= (mosaic[i].flags & ID_IMAGE_TGROUP_PHOTCAL); + image[im].flags |= (mosaic[i].flags & ID_IMAGE_MOSAIC_PHOTCAL); + image[im].flags |= (mosaic[i].flags & ID_IMAGE_IMAGE_PHOTCAL); } mosaic[i].McalPSF = 0.0; @@ -756,5 +770,5 @@ // this mosaic has been identified as poor. - mosaic[myMosaic].flags |= ID_IMAGE_PHOTOM_POOR; + mosaic[myMosaic].flags |= ID_IMAGE_MOSAIC_POOR; mosaic[myMosaic].McalPSF = 0.0; mosaic[myMosaic].McalAPER = 0.0; @@ -776,5 +790,5 @@ // this mosaic has been identified as good. - mosaic[myMosaic].flags &= ~ID_IMAGE_PHOTOM_POOR; + mosaic[myMosaic].flags &= ~ID_IMAGE_MOSAIC_POOR; image = getimages (&Nimage, NULL); @@ -1297,8 +1311,10 @@ assert (Nbright >= 0); - psfStars->alldata-> yVector[N] = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat; + float Moff = Mcal + Mgrp + Mgrid - Mflat; + + psfStars->alldata-> yVector[N] = MsysPSF - MrelPSF - Moff; psfStars->alldata->dyVector[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); - kronStars->alldata-> yVector[N] = MsysKron - MrelKron - Mcal - Mgrid + Mflat; + kronStars->alldata-> yVector[N] = MsysKron - MrelKron - Moff; kronStars->alldata->dyVector[N] = psfStars->alldata->dyVector[N]; @@ -1337,4 +1353,5 @@ myMosaic[0].McalPSF = psfStars->bSaveArray[0][0]; myMosaic[0].dMcal = psfStars->bSigma[0]; + myMosaic[0].stdev = psfStars->sigma; myMosaic[0].McalChiSq = psfStars->chisq; myMosaic[0].nFitPhotom = psfStars->Nmeas; @@ -1863,7 +1880,5 @@ void clean_mosaics () { - off_t i, N, mark, Ntotal, Nmark, Nscatter, Noffset, Nrejects, Nchisq; - double *mlist, *slist, *dlist; - double MaxOffset, MedOffset, MaxScatter; + double *mlist, *slist; if (!MOSAIC_ZEROPT) return; @@ -1874,13 +1889,12 @@ ALLOCATE (mlist, double, Nmosaic); ALLOCATE (slist, double, Nmosaic); - ALLOCATE (dlist, double, Nmosaic); - - for (i = N = 0; i < Nmosaic; i++) { + + int N = 0; + for (int i = 0; i < Nmosaic; i++) { // if (mosaic[i].flags & IMAGE_BAD) continue; if (!(mosaic[i].flags & ID_IMAGE_MOSAIC_PHOTCAL)) continue; if (mosaic[i].skipCal) continue; - mlist[N] = mosaic[i].McalPSF; - slist[N] = mosaic[i].dMcal; - dlist[N] = 1; + mlist[N] = mosaic[i].McalChiSq; + slist[N] = mosaic[i].stdev; N++; } @@ -1889,17 +1903,19 @@ liststats_setmode (&stats, "MEAN"); - liststats (mlist, dlist, NULL, N, &stats); - MaxOffset = MAX (IMAGE_OFFSET, 2*stats.sigma); - MedOffset = stats.median; - - liststats (slist, dlist, NULL, N, &stats); - MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median); - fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset); + liststats (mlist, NULL, NULL, N, &stats); + float MaxChiSq = stats.Upper90; + + // old version based on Mcal value + // MaxOffset = MAX (IMAGE_OFFSET, 2*stats.sigma); + /// MedOffset = stats.median; + + liststats (slist, NULL, NULL, N, &stats); + float MaxScatter = stats.Upper90; + // old: MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median); + + fprintf (stderr, "MOSAIC: Max ChiSq: %f, Max Scatter: %f\n", MaxChiSq, MaxScatter); - // XXX this needs to be configurable - float MaxChiSq = 10.0; - - Ntotal = Nmark = Nscatter = Noffset = Nrejects = Nchisq = 0; - for (i = 0; i < Nmosaic; i++) { + int Ntotal = 0, Nmark = 0, Nscatter = 0, Nfew = 0, Nchisq = 0; + for (int i = 0; i < Nmosaic; i++) { // if we are keeping ubercal sacrosanct, then we should not be allowed to break them... if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; @@ -1907,21 +1923,15 @@ Ntotal ++; + if (mosaic[i].skipCal) continue; + + int mark = FALSE; + if (mosaic[i].flags & (ID_IMAGE_PHOTOM_FEW)) { - continue; - } - if (mosaic[i].skipCal) continue; - - mark = FALSE; - if (mosaic[i].dMcal > MaxScatter) { + mark = TRUE; + Nfew ++; + } + if (mosaic[i].stdev > MaxScatter) { mark = TRUE; Nscatter ++; - } - if (fabs(mosaic[i].McalPSF - MedOffset) > MaxOffset) { - mark = TRUE; - Noffset ++; - } - if (mosaic[i].nFitPhotom < 0.7*mosaic[i].nValPhotom) { - mark = TRUE; - Nrejects ++; } if (mosaic[i].McalChiSq > MaxChiSq) { @@ -1937,9 +1947,8 @@ } - fprintf (stderr, OFF_T_FMT" of "OFF_T_FMT" mosaics marked poor ("OFF_T_FMT" scatter, "OFF_T_FMT" offset)\n", Nmark, Ntotal, Nscatter, Noffset); + fprintf (stderr, "%d of %d mosaics marked poor (%d scatter, %d chisq, %d few)\n", Nmark, Ntotal, Nscatter, Nchisq, Nfew); free (mlist); free (slist); - free (dlist); } Index: /trunk/Ohana/src/relphot/src/TGroupOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/TGroupOps.c (revision 41554) +++ /trunk/Ohana/src/relphot/src/TGroupOps.c (revision 41555) @@ -2,4 +2,5 @@ void assignMosaicsToTGroups (void); +int save_test_night_measures (FILE *fout, TGroup *myTGroup, Catalog *catalog, FlatCorrectionTable *flatcorr); // tgroupTimes carries the times of the tgroups (initially, just the photometric nights) @@ -82,4 +83,6 @@ tgroup[j].dMcal = 0.0; // note : at the end, tgroup.Mcal is added back to the input images tgroup[j].dMsys = 0.0; + tgroup[j].dMmin = NAN; + tgroup[j].dMmax = NAN; tgroup[j].McalChiSq = 0.0; // NAN or 0.0? @@ -142,4 +145,6 @@ tgroup[j].dMcal = 0.0; // note : at the end, tgroup.Mcal is added back to the input images tgroup[j].dMsys = 0.0; + tgroup[j].dMmin = NAN; + tgroup[j].dMmax = NAN; tgroup[j].McalChiSq = 0.0; // NAN or 0.0? @@ -421,10 +426,12 @@ off_t im = tgroup[i].image[j]; double Mgrp = (TGROUP_FIT_AIRMASS) ? tgroup[i].McalPSF + tgroup[i].dKlam*(image[im].secz - 1.0) : tgroup[i].McalPSF; - image[im].McalPSF += Mgrp; - image[im].McalAPER += Mgrp; - image[im].dMcal = tgroup[i].dMcal; - image[im].McalChiSq = tgroup[i].McalChiSq; - image[im].dMagSys = tgroup[i].dMsys; - image[im].nFitPhotom = tgroup[i].nFitPhotom; + if (tgroup[i].flags & ID_IMAGE_TGROUP_PHOTCAL) { + image[im].McalPSF = Mgrp; + image[im].McalAPER = Mgrp; + image[im].dMcal = tgroup[i].dMcal; + image[im].McalChiSq = tgroup[i].McalChiSq; + image[im].dMagSys = tgroup[i].dMsys; + image[im].nFitPhotom = tgroup[i].nFitPhotom; + } image[im].flags |= (tgroup[i].flags & ID_IMAGE_NIGHT_POOR); // probably not necessary : it is set in markBadTGroup image[im].flags |= (tgroup[i].flags & ID_IMAGE_MOSAIC_POOR); // probably not necessary : it is set in markBadTGroup @@ -432,4 +439,7 @@ image[im].flags |= (tgroup[i].flags & ID_IMAGE_PHOTOM_POOR); image[im].flags |= (tgroup[i].flags & ID_IMAGE_PHOTOM_UBERCAL); + image[im].flags |= (tgroup[i].flags & ID_IMAGE_TGROUP_PHOTCAL); + image[im].flags |= (tgroup[i].flags & ID_IMAGE_MOSAIC_PHOTCAL); + image[im].flags |= (tgroup[i].flags & ID_IMAGE_IMAGE_PHOTCAL); // fprintf (stderr, "TG to IMAGE: %f +/- %f -> %f (%d)\n", tgroup[i].McalPSF, tgroup[i].dMcal, image[im].McalPSF, tgroup[i].nFitPhotom); @@ -687,6 +697,6 @@ info->psfStars.MaxIterations = 0; } - if (1) { - FitDataSetAddPriors (&info->psfStars); + FitDataSetAddPriors (&info->psfStars); // why only psfStars? + if (TGROUP_FIT_AIRMASS) { info->psfStars.bPriorValue[1] = 0.0; // note that we are fitting relative to the nominal slope info->psfStars.bPriorSigma[1] = 0.04; // XXX this prior sigma needs to be user-configured @@ -812,6 +822,19 @@ int testImage = FALSE; TGTimes *tgroup = (TGTimes *) myTGroup->parent; - // testImage |= (abs(tgroup->start - 1245283200) < 10); + double mjdStart = ohana_sec_to_mjd (tgroup->start); + + // unset this flag at start (set below if zeropoint is calculated) + myTGroup->flags &= ~ID_IMAGE_TGROUP_PHOTCAL; + // testImage = TRUE; + testImage |= (abs(mjdStart - 57353) < 0.1); + testImage |= (abs(mjdStart - 57952) < 0.1); + testImage |= (abs(mjdStart - 55843) < 0.1); + testImage |= (abs(mjdStart - 56586) < 0.1); + testImage |= (abs(mjdStart - 57323) < 0.1); + //testImage |= (abs(mjdStart - 57971) < 0.1); + //testImage |= (abs(mjdStart - 57974) < 0.1); + //testImage |= (abs(mjdStart - 57978) < 0.1); + //testImage |= (abs(mjdStart - 57981) < 0.1); FILE *fout = NULL; @@ -820,4 +843,7 @@ snprintf (filename, 64, "test.%05d.%02d.dat", (int) Ngrp, npass_output); fout = fopen (filename, "w"); + fprintf (fout, "# tgroup %s (%.0f)\n", ohana_sec_to_date(tgroup->start), mjdStart); + save_test_night_measures (fout, myTGroup, catalog, flatcorr); + fclose (fout); } @@ -831,28 +857,4 @@ off_t c = myTGroup->catalog[j]; - // NOTE : we are only using Mcal == McalPSF in this function; we set McalAPER to McalPSF - if (fout) { - - float Mcal = getMcal (m, c, MAG_CLASS_PSF); // image zero point - float Mmos = getMmos (m, c); // mosaic zero point - // XXX these two should both be 0.0 if we are fitting tgroups : skip those calls above? - - float Mgrid = getMgrid (m, c); // camera offset (deprecated?) - float MrelPSF = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude - float Mflat = getMflat (m, c, flatcorr, catalog); // flat-field correction - - off_t n = catalog[c].measureT[m].averef; - - // magnitude for this measurement - float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); - - // for historical reasons, Mflat is defined with the wrong sign - float delta = MsysPSF - MrelPSF - Mcal - Mmos - Mgrid + Mflat; - - int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD); - - fprintf (fout, "%f %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, catalog[c].measureT[m].airmass, delta, isBad); - } - if (catalog[c].measureT[m].dbFlags & MEAS_BAD) { info->Nbad ++; @@ -916,30 +918,39 @@ assert (Nbright >= 0); - psfStars->alldata-> yVector[N] = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat; - psfStars->alldata-> xVector[N] = (catalog[c].measureT[m].airmass - 1.0); + float Moff = Mcal + Mmos + Mgrid - Mflat; + + psfStars->alldata-> yVector[N] = MsysPSF - MrelPSF - Moff; psfStars->alldata->dyVector[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); - kronStars->alldata-> yVector[N] = MsysKron - MrelKron - Mcal - Mgrid + Mflat; - kronStars->alldata-> xVector[N] = psfStars->alldata-> xVector[N]; + kronStars->alldata-> yVector[N] = MsysKron - MrelKron - Moff; kronStars->alldata->dyVector[N] = psfStars->alldata->dyVector[N]; + + if (TGROUP_FIT_AIRMASS) { + psfStars->alldata-> xVector[N] = (catalog[c].measureT[m].airmass - 1.0); + kronStars->alldata-> xVector[N] = psfStars->alldata-> xVector[N]; + } if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) { brightStars->alldata-> yVector[Nbright] = psfStars->alldata-> yVector[N]; - brightStars->alldata-> xVector[Nbright] = psfStars->alldata-> xVector[N]; brightStars->alldata->dyVector[Nbright] = psfStars->alldata->dyVector[N]; + if (TGROUP_FIT_AIRMASS) { + brightStars->alldata-> xVector[Nbright] = psfStars->alldata-> xVector[N]; + } Nbright ++; } N++; } - if (fout) { fclose (fout); } /* too few good measurements or too many bad measurements (skip in PoorImages run) */ int mark = (N < NIGHT_TOOFEW) || (N < NIGHT_GOOD_FRACTION*myTGroup->Nmeasure); if (mark) { - if (VERBOSE2) { + if (TRUE) { TGTimes *mygroup = (TGTimes *) myTGroup->parent; - fprintf (stderr, "marked tgroup %f,%d, (%d < %d) || (%d < %f*"OFF_T_FMT")\n", ohana_sec_to_mjd(mygroup->start), myTGroup->photcode, N, NIGHT_TOOFEW, N, NIGHT_GOOD_FRACTION, myTGroup->Nmeasure); + if (0) { // XXX EAM + fprintf (stderr, "marked tgroup %f,%d, (%d < %d) || (%d < %f*"OFF_T_FMT")\n", ohana_sec_to_mjd(mygroup->start), myTGroup->photcode, N, NIGHT_TOOFEW, N, NIGHT_GOOD_FRACTION, myTGroup->Nmeasure); + } } myTGroup->flags |= ID_IMAGE_PHOTOM_FEW; + myTGroup->nValPhotom = N; info->Nfew ++; if (testImage) { @@ -949,4 +960,24 @@ } else { myTGroup->flags &= ~ID_IMAGE_PHOTOM_FEW; + } + + /* we have a non-trivial fraction of measurements which are extreme outliers, and + measurements with very small errors. To avoid over-weighting measurements with very + small errors, let's add something in quadarature. + */ + + // 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) + + 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; } @@ -962,8 +993,10 @@ myTGroup->dMcal = psfStars->bSigma[0]; myTGroup->stdev = psfStars->sigma; + myTGroup->dMmin = psfStars->min; + myTGroup->dMmax = psfStars->max; myTGroup->McalChiSq = psfStars->chisq; myTGroup->nFitPhotom = psfStars->Nmeas; myTGroup->nValPhotom = N; - myTGroup->dKlam = psfStars->bSaveArray[1][0]; + myTGroup->dKlam = (TGROUP_FIT_AIRMASS) ? psfStars->bSaveArray[1][0] : 0; // XXX this does not make sense: I need to apply the airmass slope calculated above first @@ -991,4 +1024,6 @@ } + myTGroup->flags |= ID_IMAGE_TGROUP_PHOTCAL; // set this flag + if (testImage) { fprintf (stderr, "%f %f : %f\n", myTGroup->McalPSF, myTGroup->dMsys, myTGroup->McalChiSq); @@ -998,4 +1033,41 @@ } +int save_test_night_measures (FILE *fout, TGroup *myTGroup, Catalog *catalog, FlatCorrectionTable *flatcorr) { + + int Nsecfilt = GetPhotcodeNsecfilt (); + + for (int j = 0; j < myTGroup->Nmeasure; j++) { + + off_t m = myTGroup->measure[j]; + off_t c = myTGroup->catalog[j]; + + // XXX these two should both be 0.0 if we are fitting tgroups + float Mcal = getMcal (m, c, MAG_CLASS_PSF); // image zero point + float Mmos = getMmos (m, c); // mosaic zero point + + float Mgrid = getMgrid (m, c); // camera offset (deprecated?) + float MrelPSF = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude + float Mflat = getMflat (m, c, flatcorr, catalog); // flat-field correction + + off_t n = catalog[c].measureT[m].averef; + + // magnitude for this measurement + float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); + + // for historical reasons, Mflat is defined with the wrong sign + float delta = MsysPSF - MrelPSF - Mcal - Mmos - Mgrid + Mflat - myTGroup->McalPSF; + + int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD); + + double mjdMeas = ohana_sec_to_mjd (catalog[c].measureT[m].t); + + fprintf (fout, "%f %f : %f %f %f %f %f : %f %f %d : %f\n", + catalog[c].averageT[n].R, catalog[c].averageT[n].D, + MsysPSF, MrelPSF, Mcal, Mgrid, Mflat, catalog[c].measureT[m].airmass, delta, isBad, mjdMeas); + + } + return TRUE; +} + int setMgrp_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { @@ -1221,10 +1293,19 @@ } -/* mark tgroup if: abs(Mcal - ) too large, dMcal too large */ +/* new rules for bad tgroups (bad nights) + * too few exposures (Nmosaic < X) + * dMsys or stdev > X + * X%-ile of dMsys or stdev + * McalChiSq > X + * X%-ile of McalChiSq + + * I can define rules for these but they would be arbitrary + * I think I should do an analysis of the observed distribution + + */ + void clean_tgroups () { - off_t mark, Ntotal, Nmark, Nscatter, Noffset, Nrejects, Nchisq; - double *mlist, *slist, *dlist; - double MaxOffset, MedOffset, MaxScatter; + double *mlist, *slist; if (!TGROUP_ZEROPT) return; @@ -1234,5 +1315,4 @@ ALLOCATE (mlist, double, NtgroupTimes*Nphotcodes); ALLOCATE (slist, double, NtgroupTimes*Nphotcodes); - ALLOCATE (dlist, double, NtgroupTimes*Nphotcodes); int N = 0; @@ -1241,7 +1321,7 @@ for (int j = 0; j < tgroupTimes[i][0].nCode; j++) { if (tgroup[j].flags & IMAGE_BAD) continue; - mlist[N] = tgroup[j].McalPSF; - slist[N] = tgroup[j].dMcal; - dlist[N] = 1; + mlist[N] = tgroup[j].McalChiSq; + slist[N] = tgroup[j].stdev; // stdev of all measurements + // slist[N] = tgroup[j].dMsys; // stdev of bright measurements N++; } @@ -1251,16 +1331,13 @@ liststats_setmode (&stats, "MEAN"); - liststats (mlist, dlist, NULL, N, &stats); - MaxOffset = MAX (NIGHT_OFFSET, 3*stats.sigma); - MedOffset = stats.median; - - liststats (slist, dlist, NULL, N, &stats); - MaxScatter = MAX (NIGHT_SCATTER, 2*stats.median); - fprintf (stderr, "TGROUPS: median Mrel: %f, median dMrel: %f, Max Offset: %f, Max Scatter: %f\n", MedOffset, stats.median, MaxOffset, MaxScatter); + liststats (mlist, NULL, NULL, N, &stats); + float MaxChiSq = stats.Upper90; + + liststats (slist, NULL, NULL, N, &stats); + float MaxScatter = stats.Upper90; + + fprintf (stderr, "TGROUPS: Max ChiSq: %f, Max Scatter: %f\n", MaxChiSq, MaxScatter); - // XXX this needs to be configurable - float MaxChiSq = 10.0; - - Ntotal = Nmark = Nscatter = Noffset = Nrejects = Nchisq = 0; + int Ntotal = 0, Nmark = 0, Nscatter = 0, Nchisq = 0, NfewNights = 0, NfewExp = 0; for (int i = 0; i < NtgroupTimes; i++) { TGroup *tgroup = tgroupTimes[i][0].byCode; @@ -1269,25 +1346,22 @@ Ntotal ++; - // ignore the nights with too few measurements - // XXX what to do with nights with too few measurements? - // XXX probably free the images to be fitted + int mark = FALSE; + + // too few measurements (set in setMgrp_tgroup) if (tgroup[j].flags & ID_IMAGE_PHOTOM_FEW) { - markBadTGroup(&tgroup[j]); // mark the image associated with a bad night - continue; + mark = TRUE; + NfewNights ++; } - - mark = FALSE; - if (tgroup[j].dMcal > MaxScatter) { + // too few exposures (configure) + if (tgroup[j].Nmosaic < 4) { + mark = TRUE; + NfewExp ++; + } + // scatter too large + if (tgroup[j].stdev > MaxScatter) { mark = TRUE; Nscatter ++; } - if (fabs(tgroup[j].McalPSF - MedOffset) > MaxOffset) { - mark = TRUE; - Noffset ++; - } - if (tgroup[j].nFitPhotom < 0.7*tgroup[j].nValPhotom) { - mark = TRUE; - Nrejects ++; - } + // chisq too large if (tgroup[j].McalChiSq > MaxChiSq) { mark = TRUE; @@ -1303,9 +1377,8 @@ } - fprintf (stderr, OFF_T_FMT" of "OFF_T_FMT" tgroups marked poor ("OFF_T_FMT" scatter, "OFF_T_FMT" offset, "OFF_T_FMT" rejects, "OFF_T_FMT" chisq)\n", Nmark, Ntotal, Nscatter, Noffset, Nrejects, Nchisq); + fprintf (stderr, "%d of %d tgroups marked poor (%d scatter, %d few nights, %d few exposures, %d chisq)\n", Nmark, Ntotal, Nscatter, NfewNights, NfewExp, Nchisq); free (mlist); free (slist); - free (dlist); } @@ -1465,2 +1538,70 @@ } } + +void dump_tgroup_imstats (int Npass) { + + FILE *fout = NULL; + char filename[64]; + + snprintf (filename, 64, "imstats.grp.%02d.dat", Npass); + fout = fopen (filename, "w"); + + for (off_t i = 0; i < NtgroupTimes; i++) { + + TGroup *tgroup = tgroupTimes[i][0].byCode; + for (off_t j = 0; j < tgroupTimes[i][0].nCode; j++) { + fprintf (fout, "%.0f %5d | %6ld %6d %6d %4ld %4ld | 0x%08x %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\n", + ohana_sec_to_mjd(tgroupTimes[i][0].start), tgroup[j].photcode, + tgroup[j].Nmeasure, tgroup[j].nFitPhotom, tgroup[j].nValPhotom, tgroup[j].Nimage, tgroup[j].Nmosaic, + tgroup[j].flags, tgroup[j].McalPSF, tgroup[j].dKlam, tgroup[j].dMcal, tgroup[j].dMsys, tgroup[j].stdev, tgroup[j].dMmin, tgroup[j].dMmax, tgroup[j].McalChiSq); + } + } + fclose (fout); + + // ************ mosaic stats + + snprintf (filename, 64, "imstats.mos.%02d.dat", Npass); + fout = fopen (filename, "w"); + + for (off_t i = 0; i < NtgroupTimes; i++) { + + TGroup *tgroup = tgroupTimes[i][0].byCode; + for (off_t j = 0; j < tgroupTimes[i][0].nCode; j++) { + + for (off_t k = 0; k < tgroup->Nmosaic; k++) { + Mosaic *mosaic = tgroup[j].mosaic[k]; + + fprintf (fout, "%4ld %3ld %3ld %.0f %5d %5d %5d | %12.5f %5.3f 0x%08x | %6.3f %6.3f %6.3f %6.3f\n", + i, j, k, ohana_sec_to_mjd(tgroupTimes[i][0].start), mosaic->photcode, mosaic->nFitPhotom, mosaic->nValPhotom, + ohana_sec_to_mjd(mosaic->start), mosaic->secz, mosaic->flags, + mosaic->McalPSF, mosaic->dMcal, mosaic->dMsys, mosaic->McalChiSq); + } + } + } + fclose (fout); + + // ************ image stats + + snprintf (filename, 64, "imstats.chp.%02d.dat", Npass); + fout = fopen (filename, "w"); + + off_t Nimage = 0; + Image *image = getimages (&Nimage, NULL); + for (off_t i = 0; i < NtgroupTimes; i++) { + + TGroup *tgroup = tgroupTimes[i][0].byCode; + for (off_t j = 0; j < tgroupTimes[i][0].nCode; j++) { + + for (off_t k = 0; k < tgroup->Nimage; k++) { + off_t im = tgroup[j].image[k]; + + fprintf (fout, "%4ld %3ld %3ld %.0f %5d | %12.5f %5.3f 0x%08x | %6.3f %6.3f %6.3f %6.3f %6.3f\n", + i, j, k, ohana_sec_to_mjd(tgroupTimes[i][0].start), image[im].photcode, + ohana_sec_to_mjd(image[im].tzero), image[im].secz, image[im].flags, + image[im].McalPSF, image[im].McalAPER, image[im].dMcal, image[im].dMagSys, image[im].McalChiSq); + } + } + } + + fclose (fout); +} Index: /trunk/Ohana/src/relphot/src/ZeroPointModes.c =================================================================== --- /trunk/Ohana/src/relphot/src/ZeroPointModes.c (revision 41554) +++ /trunk/Ohana/src/relphot/src/ZeroPointModes.c (revision 41555) @@ -4,4 +4,15 @@ void SetZptIteration (int current) { + // XXX add 'manual' option here: ask for desired zpt iteration + if (MANUAL_ITERATION) { + int newValue; + fprintf (stdout, "Set Iteration (current = %d): ", CurrentLoop); + fscanf (stdin, "%d", &newValue); + fprintf (stdout, "Using iteration %d\n", newValue); + CurrentLoop = newValue; + return; + } + + // automatic iteration: CurrentLoop = current; } @@ -14,5 +25,10 @@ if (TGROUP_ZEROPT && MOSAIC_ZEROPT) { - if ((CurrentLoop > 4) && (CurrentLoop <= 8)) { + if (CurrentLoop < 6) { + TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_ALL; + IMAGE_ZPT_MODE = IMAGE_ZPT_MODE_NONE; + MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE; + } + if ((CurrentLoop >= 6) && (CurrentLoop < 12)) { // after iterating a few times on the TGroup zero points: // * identify the photometric nights @@ -23,5 +39,5 @@ clean_tgroups(); // do this on every pass to update the status of nights } - if ((CurrentLoop > 8) && (CurrentLoop <= 999)) { + if ((CurrentLoop >= 12) && (CurrentLoop < 18)) { // after iterating a few times on the TGroup & Mosaic zero points: // * identify the good / bad Mosaics @@ -33,7 +49,23 @@ clean_mosaics(); // do this on every pass } + if ((CurrentLoop >= 18) && (CurrentLoop <= 999)) { + // after iterating a few times on the TGroup & Mosaic zero points: + // * identify the good / bad Mosaics + // * fit the chips from the bad Mosaics + TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_GOOD_NIGHT; // stop fitting the bad nights + MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_BAD_NIGHT_GOOD_MOSAIC; // only fit good mosaics on bad nights + IMAGE_ZPT_MODE = IMAGE_ZPT_MODE_BAD_NIGHT_BAD_MOSAIC; // only fit bad mosaics on bad nights + clean_tgroups(); // do this on every pass + clean_mosaics(); // do this on every pass + clean_images(); // do this on every pass + } } if (TGROUP_ZEROPT && !MOSAIC_ZEROPT) { + if (CurrentLoop <= 4) { + TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_ALL; + IMAGE_ZPT_MODE = IMAGE_ZPT_MODE_NONE; + MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE; + } if ((CurrentLoop > 4) && (CurrentLoop <= 999)) { // after iterating a few times on the TGroup zero points: @@ -48,4 +80,9 @@ if (!TGROUP_ZEROPT && MOSAIC_ZEROPT) { + if (CurrentLoop <= 4) { + TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_NONE; + MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_ALL; + IMAGE_ZPT_MODE = IMAGE_ZPT_MODE_NONE; + } if ((CurrentLoop > 4) && (CurrentLoop <= 999)) { // after iterating a few times on the Mosaic zero points: @@ -60,11 +97,8 @@ if (!TGROUP_ZEROPT && !MOSAIC_ZEROPT) { - if ((CurrentLoop > 4) && (CurrentLoop <= 999)) { - // after iterating a few times on the Image zero points: - // * ??? - TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_NONE; - MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE; - IMAGE_ZPT_MODE = IMAGE_ZPT_MODE_ALL; - } + // if we are not fitting TGROUP or MOSAIC, just fit all images + TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_NONE; + MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE; + IMAGE_ZPT_MODE = IMAGE_ZPT_MODE_ALL; } } @@ -83,14 +117,30 @@ return FALSE; case ZPT_MOSAIC: + if (CurrentLoop < 9) return TRUE; + return FALSE; + case ZPT_IMAGES: + default: + if (CurrentLoop < 15) return TRUE; + return FALSE; + } + } + + if (TGROUP_ZEROPT && !MOSAIC_ZEROPT) { + switch (mode) { + case ZPT_STARS: + if (CurrentLoop < 3) return TRUE; + return FALSE; + case ZPT_TGROUP: + case ZPT_MOSAIC: + if (CurrentLoop < 3) return TRUE; + return FALSE; + case ZPT_IMAGES: + default: if (CurrentLoop < 7) return TRUE; return FALSE; - case ZPT_IMAGES: - default: - if (CurrentLoop < 11) return TRUE; - return FALSE; - } - } - - if (TGROUP_ZEROPT && !MOSAIC_ZEROPT) { + } + } + + if (!TGROUP_ZEROPT && MOSAIC_ZEROPT) { switch (mode) { case ZPT_STARS: @@ -108,20 +158,4 @@ } - if (!TGROUP_ZEROPT && MOSAIC_ZEROPT) { - switch (mode) { - case ZPT_STARS: - if (CurrentLoop < 3) return TRUE; - return FALSE; - case ZPT_TGROUP: - case ZPT_MOSAIC: - if (CurrentLoop < 3) return TRUE; - return FALSE; - case ZPT_IMAGES: - default: - if (CurrentLoop < 7) return TRUE; - return FALSE; - } - } - if (!TGROUP_ZEROPT && !MOSAIC_ZEROPT) { switch (mode) { @@ -139,2 +173,30 @@ return FALSE; } + +/* TGROUP & MOSAIC: + 0 - TGROUP ONLY & OLS + 1 - TGROUP ONLY & OLS + 2 - TGROUP ONLY & OLS + + 3 - TGROUP ONLY & IRLS + 4 - TGROUP ONLY & IRLS + 5 - TGROUP ONLY & IRLS + + 6 - TGROUP + MOSAIC & OLS MOSAIC + 7 - TGROUP + MOSAIC & OLS MOSAIC + 8 - TGROUP + MOSAIC & OLS MOSAIC + + 9 - TGROUP + MOSAIC & IRLS MOSAIC + 10 - TGROUP + MOSAIC & IRLS MOSAIC + 11 - TGROUP + MOSAIC & IRLS MOSAIC + + 12 - TGROUP + MOSAIC + IMAGE & OLS IMAGE + 13 - TGROUP + MOSAIC + IMAGE & OLS IMAGE + 14 - TGROUP + MOSAIC + IMAGE & OLS IMAGE + + 15 - TGROUP + MOSAIC + IMAGE & IRLS IMAGE + 16 - TGROUP + MOSAIC + IMAGE & IRLS IMAGE + 17 - TGROUP + MOSAIC + IMAGE & IRLS IMAGE + +*/ + Index: /trunk/Ohana/src/relphot/src/args.c =================================================================== --- /trunk/Ohana/src/relphot/src/args.c (revision 41554) +++ /trunk/Ohana/src/relphot/src/args.c (revision 41555) @@ -192,4 +192,10 @@ } + MANUAL_ITERATION = FALSE; + if ((N = get_argument (argc, argv, "-manual-iteration"))) { + remove_argument (N, &argc, argv); + MANUAL_ITERATION = TRUE; + } + NGRID = 8; if ((N = get_argument (argc, argv, "-ngrid"))) { @@ -325,5 +331,4 @@ USE_MCAL_PSF_FOR_STACK_APER = TRUE; } - // default behavior is to fit all chips independently Index: /trunk/Ohana/src/relphot/src/liststats.c =================================================================== --- /trunk/Ohana/src/relphot/src/liststats.c (revision 41554) +++ /trunk/Ohana/src/relphot/src/liststats.c (revision 41555) @@ -26,4 +26,6 @@ stats->Upper80 = NAN; stats->Lower20 = NAN; + stats->Upper90 = NAN; + stats->Lower10 = NAN; stats->total = NAN; stats->Nmeas = 0; @@ -49,6 +51,8 @@ if (N < 1) return (FALSE); + int N90 = MIN (N-1, 0.9*N); int N80 = MIN (N-1, 0.8*N); int N20 = MAX (0, 0.2*N); + int N10 = MAX (0, 0.1*N); if (weight) { @@ -68,4 +72,6 @@ stats[0].Upper80 = value[N80]; stats[0].Lower20 = value[N20]; + stats[0].Upper90 = value[N90]; + stats[0].Lower10 = value[N10]; switch (stats->statmode) { @@ -116,4 +122,5 @@ WeightedMean |= (stats->statmode == STATS_INNER_WTMEAN); WeightedMean |= (stats->statmode == STATS_CHI_INNER_WTMEAN); + if (!dvalue) WeightedMean = FALSE; // warn the user? /* calculating the per-star offset based on the desired weighting scheme */ @@ -137,10 +144,10 @@ M += value[i] * weight[i]; W += weight[i]; - dM += SQ (weight[i] * dvalue[i]); + dM += dvalue ? SQ (weight[i] * dvalue[i]) : 0.0; R += weight[i]; Nm += 1.0; } Mo = M / W; - dMo = sqrt (dM) / R; + dMo = dvalue ? sqrt (dM) / R : NAN; } } else { @@ -159,9 +166,9 @@ for (i = ks; i < ke; i++) { M += value[i]; - dM += SQ (dvalue[i]); + dM += dvalue ? SQ (dvalue[i]) : 0.0; Nm += 1.0; } Mo = M / Nm; - dMo = sqrt (dM) / Nm; + dMo = dvalue ? sqrt (dM) / Nm : NAN; } } Index: /trunk/Ohana/src/relphot/src/load_images.c =================================================================== --- /trunk/Ohana/src/relphot/src/load_images.c (revision 41554) +++ /trunk/Ohana/src/relphot/src/load_images.c (revision 41555) @@ -58,7 +58,6 @@ // 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); - } + // reset these two regardless (we will re-determine) + subset[i].flags &= ~(ID_IMAGE_PHOTOM_NOCAL | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP | ID_IMAGE_PHOTOM_FEW); if (RESET) { if (RESET_ZEROPTS) { @@ -70,5 +69,6 @@ } } - subset[i].flags &= ~ID_IMAGE_PHOTOM_POOR; + subset[i].flags &= ~(ID_IMAGE_MOSAIC_POOR | ID_IMAGE_NIGHT_POOR); + subset[i].flags &= ~(ID_IMAGE_IMAGE_PHOTCAL | ID_IMAGE_MOSAIC_PHOTCAL | ID_IMAGE_TGROUP_PHOTCAL); subset[i].ubercalDist = 1000; } Index: /trunk/Ohana/src/relphot/src/relphot_images.c =================================================================== --- /trunk/Ohana/src/relphot/src/relphot_images.c (revision 41554) +++ /trunk/Ohana/src/relphot/src/relphot_images.c (revision 41555) @@ -117,10 +117,15 @@ for (i = 0; i < NLOOP; i++) { SetZptIteration (i); + dump_tgroups (catalog, i); dump_catalog (catalog, 0, i); // for a test, just dump a specific catalog + dump_tgroup_imstats (i); setMrel (catalog, Ncatalog, flatcorr); // threaded + setMcalTest (catalog, flatcorr); + dump_tgroups (catalog, i + 20); dump_catalog (catalog, 0, i + 20); // for a test, just dump a specific catalog + dump_tgroup_imstats (i + 20); setMcal (catalog, FALSE, flatcorr); @@ -129,4 +134,5 @@ dump_tgroups (catalog, i + 40); dump_catalog (catalog, 0, i + 40); // for a test, just dump a specific catalog + dump_tgroup_imstats (i + 40); setMgrid (catalog, flatcorr); @@ -143,6 +149,6 @@ global_stats (catalog, Ncatalog, flatcorr, i); - SetZeroPointModes (); + MARKTIME("-- finished loop %d: %f sec\n", i, dtime); fprintf (stderr, "----- continue loop %d -----\n", i); Index: /trunk/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- /trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 41554) +++ /trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 41555) @@ -272,4 +272,8 @@ float Map = NAN; float dMap = NAN; + + // XXX I am not convinced the sign on Mflat here is correct + float Moff = Mflat + Mcal + Mmos + Mgrp + Mgrid; + if (isSetMrelFinal) { Map = PhotCat (&measure[k], MAG_CLASS_APER); @@ -282,5 +286,5 @@ if (CHECK_VALID_MAG(Map, dMap)) { int Nap = results->aperData[Nsec].Nlist; - results->aperData[Nsec].flxlist[Nap] = Map - Mflat - Mcal - Mmos - Mgrp - Mgrid; // this is consistent with PhotRel + results->aperData[Nsec].flxlist[Nap] = Map - Moff; // this is consistent with PhotRel results->aperData[Nsec].errlist[Nap] = dMap; results->aperData[Nsec].wgtlist[Nap] = modifiedWeight; @@ -302,5 +306,5 @@ if (CHECK_VALID_MAG(Mkron, dMkron)) { int Nkron = results->kronData[Nsec].Nlist; - results->kronData[Nsec].flxlist[Nkron] = Mkron - Mflat - Mcal - Mmos - Mgrp - Mgrid; // this is consistent with PhotRel + results->kronData[Nsec].flxlist[Nkron] = Mkron - Moff; // this is consistent with PhotRel results->kronData[Nsec].errlist[Nkron] = dMkron; results->kronData[Nsec].wgtlist[Nkron] = modifiedWeight; @@ -315,5 +319,5 @@ if (CHECK_VALID_MAG(Mpsf, dMpsf)) { int Npsf = results->psfData[Nsec].Nlist; - results->psfData[Nsec].flxlist[Npsf] = Mpsf - Mflat - Mcal - Mmos - Mgrp - Mgrid; // this is consistent with PhotRel + results->psfData[Nsec].flxlist[Npsf] = Mpsf - Moff; // this is consistent with PhotRel results->psfData[Nsec].errlist[Npsf] = dMpsf; results->psfData[Nsec].wgtlist[Npsf] = modifiedWeight;