Index: /trunk/Ohana/src/relphot/Makefile =================================================================== --- /trunk/Ohana/src/relphot/Makefile (revision 41555) +++ /trunk/Ohana/src/relphot/Makefile (revision 41556) @@ -54,4 +54,6 @@ $(SRC)/load_catalogs.$(ARCH).o \ $(SRC)/load_images.$(ARCH).o \ +$(SRC)/save_images.$(ARCH).o \ +$(SRC)/MagResidSave.$(ARCH).o \ $(SRC)/synthetic_mags.$(ARCH).o \ $(SRC)/plot_scatter.$(ARCH).o \ @@ -117,4 +119,6 @@ $(SRC)/select_images.$(ARCH).o \ $(SRC)/load_images.$(ARCH).o \ +$(SRC)/save_images.$(ARCH).o \ +$(SRC)/MagResidSave.$(ARCH).o \ $(SRC)/load_catalogs.$(ARCH).o \ $(SRC)/reload_catalogs.$(ARCH).o \ Index: /trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- /trunk/Ohana/src/relphot/include/relphot.h (revision 41555) +++ /trunk/Ohana/src/relphot/include/relphot.h (revision 41556) @@ -643,7 +643,7 @@ int setMcalOutput PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); -void setMcal PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr)); -int setMmos PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr)); -int setMgrp PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr)); +void setMcal PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); +int setMmos PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); +int setMgrp PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); int setMrel PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); @@ -841,6 +841,4 @@ void RepairWarpMeasuresCloseLogfile (); -int save_images_backup (FITS_DB *db); - void dump_tgroups (Catalog *catalog, int Npass); void dump_catalog (Catalog *catalog, off_t c, int Npass); @@ -853,2 +851,7 @@ int UseStandardOLS (ZptFitModeType mode); int GetActivePhotcodeIndex (int photcode); + +int save_images_backup (FITS_DB *db); +int save_images_updates (FITS_DB *db); + +int MagResidSave(char *filename, Catalog *catalog); Index: /trunk/Ohana/src/relphot/src/ImageOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/ImageOps.c (revision 41555) +++ /trunk/Ohana/src/relphot/src/ImageOps.c (revision 41556) @@ -75,4 +75,16 @@ Image *getimage (off_t N) { return (&image[N]); +} + +off_t *get_N_onImage() { + return N_onImage; +} + +IDX_T **get_ImageToCatalog() { + return ImageToCatalog; +} + +IDX_T **get_ImageToMeasure() { + return ImageToMeasure; } @@ -570,8 +582,8 @@ /* determine Mcal values for all images */ -void setMcal (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { +void setMcal (Catalog *catalog, FlatCorrectionTable *flatcorr) { off_t i, j, m, c, n; - int mark, bad, Nfew, Nbad, Nmos, Ngrp, Nrel, Ngrid, Nsys; + int mark, Nfew, Nbad, Nmos, Ngrp, Nrel, Ngrid, Nsys; if (IMAGE_ZPT_MODE == IMAGE_ZPT_MODE_NONE) return; @@ -580,10 +592,4 @@ int Nsecfilt = GetPhotcodeNsecfilt (); - - if (PoorImages) { - // XXX use bad stars and measurements for PoorImages? or not? - // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0; - IMAGE_BAD = 0; - } off_t Nmax = 0; @@ -618,10 +624,4 @@ } - /* on PoorImages run, skip good images */ - if (PoorImages) { - bad = image[i].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP); - if (!bad) continue; - } - if (image[i].photcode == 0) continue; // skip the PHU images @@ -629,4 +629,5 @@ int badMosaic = (image[i].flags & ID_IMAGE_MOSAIC_POOR); int isMosaic = isMosaicChip(image[i].photcode); + int useMcal = TRUE; // unset at start by default (set if actually used below) @@ -638,16 +639,16 @@ // in BAD_NIGHT mode, we fit ONLY images in bad nights if ((IMAGE_ZPT_MODE == IMAGE_ZPT_MODE_BAD_NIGHT)) { - if (!badNight) continue; + if (!badNight) useMcal = FALSE; // do not fit images from good nights } // in BAD_NIGHT_BAD_MOSAIC mode, we fit ONLY bad mosaics in bad nights if ((IMAGE_ZPT_MODE == IMAGE_ZPT_MODE_BAD_NIGHT_BAD_MOSAIC)) { - if (!badNight) continue; - if (!badMosaic && isMosaic) continue; // skip good mosaics (but not good images) + if (!badNight) useMcal = FALSE; // do not fit images from good nights + if (!badMosaic && isMosaic) useMcal = FALSE; // skip good mosaics (but not good images) } // in BAD_MOSAIC mode, we fit bad mosaics ignoring night state if ((IMAGE_ZPT_MODE == IMAGE_ZPT_MODE_BAD_MOSAIC)) { - if (!badMosaic && isMosaic) continue; // skip good mosaics (but not bad images) + if (!badMosaic && isMosaic) useMcal = FALSE; // skip good mosaics (but not bad images) } @@ -759,19 +760,48 @@ /* too few good measurements or too many bad measurements */ - mark = FALSE; - if (!PoorImages) { - mark = (Nref < IMAGE_TOOFEW) || (Nref < IMAGE_GOOD_FRACTION*N_onImage[i]); - if (mark) { - image[i].flags |= ID_IMAGE_PHOTOM_FEW; - Nfew ++; - } else { - image[i].flags &= ~ID_IMAGE_PHOTOM_FEW; - } - } + mark = (Nref < IMAGE_TOOFEW) || (Nref < IMAGE_GOOD_FRACTION*N_onImage[i]); + if (mark) { + image[i].flags |= ID_IMAGE_PHOTOM_FEW; + image[i].dMcal = NAN; + image[i].McalChiSq = NAN; + image[i].nFitPhotom = 0; + Nfew ++; + } else { + image[i].flags &= ~ID_IMAGE_PHOTOM_FEW; + } if (mark) continue; + // use liststats to find the 20-pct, median, 80-pct points + StatType stats; + liststats_setmode (&stats, "MEDIAN"); + liststats (psfStars.alldata->yVector, NULL, NULL, Nref, &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 < Nref; j++) { + double newSigma = hypot(psfStars.alldata->dyVector[j], 0.1*altSigma); + psfStars.alldata->dyVector[j] = newSigma; + } + // 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? + fit1d_irls (&psfStars, Nref); - image[i].McalPSF = psfStars.bSaveArray[0][0]; + fit1d_irls (&kronStars, Nkron); + fit1d_irls (&brightStars, Nbright); + + // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not) + if (useMcal) { + image[i].McalPSF = psfStars.bSaveArray[0][0]; + image[i].McalAPER = kronStars.bSaveArray[0][0]; + image[i].flags |= ID_IMAGE_IMAGE_PHOTCAL; // set this flag (unset by default) + } else { + image[i].McalPSF = 0.0; + image[i].McalAPER = 0.0; + } + + // record statistics on the fit regardless if it is used or not image[i].dMcal = psfStars.bSigma[0]; image[i].nFitPhotom = psfStars.Nmeas; @@ -779,11 +809,6 @@ Ncalibrated ++; - // 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 EAM : temporarily use this to track calibration quality - if (0) { - fit1d_irls (&kronStars, Nkron); - image[i].McalAPER = kronStars.bSaveArray[0][0]; - } + // bright end scatter + image[i].dMagSys = brightStars.sigma; if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) { @@ -796,14 +821,4 @@ } - // bright end scatter - fit1d_irls (&brightStars, Nbright); - image[i].dMagSys = brightStars.sigma; - - if (image[i].McalPSF < -CLOUD_TOLERANCE) { - 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 // since this particular image is tied to that star, bump its distance by 1 @@ -813,10 +828,4 @@ 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); - - if (PoorImages) { - IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP; - STAR_BAD = ID_OBJ_POOR | ID_OBJ_FEW; - MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA; - } FitDataSetFree (&brightStars); @@ -1291,17 +1300,2 @@ 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/MagResidSave.c =================================================================== --- /trunk/Ohana/src/relphot/src/MagResidSave.c (revision 41556) +++ /trunk/Ohana/src/relphot/src/MagResidSave.c (revision 41556) @@ -0,0 +1,253 @@ +# include "relphot.h" + +// STATUS is value expected for success +# define CHECK_STATUS(STATUS,MSG,...) \ + if (!(STATUS)) { \ + fprintf (stderr, MSG, __VA_ARGS__); \ + return FALSE; \ + } + +typedef struct { + float Mresid; + float dMresid; + float Mchisq; + int Nfit; +} MagResidResult; + +MagResidResult MagResidCalc (Catalog *catalog, Image *image, FitDataSet *psfStars, int Nsecfilt, int N_onImage, IDX_T *ImageToMeasure, IDX_T *ImageToCatalog); + +off_t *get_N_onImage(); +IDX_T **get_ImageToCatalog(); +IDX_T **get_ImageToMeasure(); + +int MagResidSave(char *filename, Catalog *catalog) { + + int i; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + off_t Nimage; + Image *image = getimages (&Nimage, NULL); + off_t *N_onImage = get_N_onImage(); + + // find the maximum number of stars (needed to allocate the psfStars structure) + 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, using OLS + FitDataSet psfStars; + FitDataSetAlloc (&psfStars, Nmax, 0, 0); + psfStars.MaxIterations = 0; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + gfits_create_table_header (&theader, "BINTABLE", "MAG_PSF_RESID"); + + gfits_define_bintable_column (&theader, "D", "RA", "ra", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "DEC", "dec", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_RESID", "magnitude residual", "magnitudes", 1.0, 0.0); + 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", "AIRMASS", "", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "MJD", "", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "EXPTIME", "", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "PHOTCODE", "", "unitless", 1.0, 0.0); + 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); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // create intermediate storage arrays + ALLOCATE_PTR (Rs, double, Nimage); + ALLOCATE_PTR (Ds, double, Nimage); + + ALLOCATE_PTR (Mresid, float, Nimage); + ALLOCATE_PTR (dMresid, float, Nimage); + ALLOCATE_PTR (Mchisq, float, Nimage); + ALLOCATE_PTR (airmass, float, Nimage); + ALLOCATE_PTR (mjd, double, Nimage); + ALLOCATE_PTR (exptime, float, Nimage); + ALLOCATE_PTR (photcode, int, Nimage); + ALLOCATE_PTR (flags, int, Nimage); + ALLOCATE_PTR (Nfit, int, Nimage); + + // XXXX fill in the vectors here + int Nsecfilt = GetPhotcodeNsecfilt (); + IDX_T **ImageToCatalog = get_ImageToCatalog(); + IDX_T **ImageToMeasure = get_ImageToMeasure(); + + for (i = 0; i < Nimage; i++) { + double Ro, Do; + XY_to_RD (&Ro, &Do, (double) 0.5*image[i].NX, (double) 0.5*image[i].NY, &image[i].coords); + + Rs[i] = Ro; + Ds[i] = Do; + Mresid[i] = NAN; + Mchisq[i] = NAN; + dMresid[i] = NAN; + Nfit[i] = 0; + Nmeas[i] = N_onImage[i]; + airmass[i] = image[i].secz; + mjd[i] = ohana_sec_to_mjd (image[i].tzero); + exptime[i] = image[i].exptime; + photcode[i] = image[i].photcode; + flags[i] = image[i].flags; + + if (image[i].photcode == 0) continue; // skip the PHU images + if (N_onImage[i] == 0) continue; + + MagResidResult result = MagResidCalc (catalog, &image[i], &psfStars, Nsecfilt, N_onImage[i], ImageToMeasure[i], ImageToCatalog[i]); + + Mresid[i] = result.Mresid; + Mchisq[i] = result.Mchisq; + dMresid[i] = result.dMresid; + 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, "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, "AIRMASS", airmass, Nimage); + gfits_set_bintable_column (&theader, &ftable, "MJD", mjd, Nimage); + gfits_set_bintable_column (&theader, &ftable, "EXPTIME", exptime, Nimage); + gfits_set_bintable_column (&theader, &ftable, "PHOTCODE", photcode, Nimage); + gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags, Nimage); + gfits_set_bintable_column (&theader, &ftable, "NFIT", Nfit, Nimage); + + free (Rs); + free (Ds); + free (Mresid ); + free (dMresid ); + free (Mchisq ); + free (airmass ); + free (mjd ); + free (exptime ); + free (photcode); + free (flags ); + free (Nfit ); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open meanmag file for output %s\n", filename); + return FALSE; + } + + int status; + status = gfits_fwrite_header (f, &header); + CHECK_STATUS (status, "ERROR: cannot write header for meanmags %s\n", filename); + + status = gfits_fwrite_matrix (f, &matrix); + CHECK_STATUS (status, "ERROR: cannot write matrix for meanmags %s\n", filename); + + status = gfits_fwrite_Theader (f, &theader); + CHECK_STATUS (status, "ERROR: cannot write table header for meanmags %s\n", filename); + + status = gfits_fwrite_table (f, &ftable); + CHECK_STATUS (status, "ERROR: cannot write table data for meanmags %s\n", filename); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + int fd = fileno (f); + + status = fflush (f); + CHECK_STATUS (!status, "ERROR: cannot flush file meanmags %s\n", filename); + + status = fsync (fd); + CHECK_STATUS (!status, "ERROR: cannot flush file meanmags %s\n", filename); + + status = fclose (f); + CHECK_STATUS (!status, "ERROR: problem closing meanmags file %s\n", filename); + + return TRUE; +} + +// determine Mag residual for this image +MagResidResult MagResidCalc (Catalog *catalog, Image *image, FitDataSet *psfStars, int Nsecfilt, int N_onImage, IDX_T *ImageToMeasure, IDX_T *ImageToCatalog) { + + off_t Nref = 0; // number of stars used to measure McalPSF + + MagResidResult result; + result.Mresid = NAN; + result.dMresid = NAN; + result.Mchisq = NAN; + result.Nfit = 0; + + for (off_t j = 0; j < N_onImage; j++) { + + off_t m = ImageToMeasure[j]; + off_t c = ImageToCatalog[j]; + + // XXX drop this? + if (catalog[c].measureT[m].dbFlags & MEAS_BAD) continue; + + float Mcal = getMcal (m, c, MAG_CLASS_PSF); + if (isnan(Mcal)) continue; + + float Mmos = getMmos (m, c); + if (isnan(Mmos)) continue; + + float Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass, NULL); // ignore error for now? + if (isnan(Mgrp)) continue; + + float Mgrid = getMgrid (m, c); + if (isnan(Mgrid)) 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); + + // MrelPSF is the average magnitude for this star. + float MrelPSF = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); + if (isnan(MrelPSF)) continue; + + // get the PSF magnitude for thie measurement, with airmass slope applied + off_t n = catalog[c].measureT[m].averef; + + // MsysPSF is the specific measurement for this star on this image + float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); + if (isnan(MsysPSF)) continue; + + // XXX skip Mflat for now + // float Moff = Mcal + Mmos + Mgrp + Mgrid - Mflat; + + float Moff = Mcal + Mmos + Mgrp + Mgrid; + + psfStars->alldata-> yVector[Nref] = MsysPSF - MrelPSF - Moff; + psfStars->alldata->dyVector[Nref] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); + Nref++; + } + + if (Nref < 5) { + result.Mresid = NAN; + result.dMresid = NAN; + result.Mchisq = NAN; + result.Nfit = Nref; + return result; + } + + // 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); + result.Mresid = psfStars->bSaveArray[0][0]; + result.dMresid = psfStars->bSigma[0]; + result.Nfit = psfStars->Nmeas; + result.Mchisq = psfStars->chisq; + + return result; +} Index: /trunk/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/MosaicOps.c (revision 41555) +++ /trunk/Ohana/src/relphot/src/MosaicOps.c (revision 41556) @@ -1,4 +1,5 @@ # include "relphot.h" off_t findMosaic (unsigned int *startTimes, off_t Nmosaic, unsigned int start); +int save_test_mosaic_measures (FILE *fout, int Nmos, Catalog *catalog, FlatCorrectionTable *flatcorr); // see discussion in ImagesOps.c re: IDX_T @@ -983,5 +984,4 @@ int Nskip; off_t Nmax; - int PoorImages; FitDataSet psfStars; FitDataSet kronStars; @@ -1002,7 +1002,7 @@ int setMmos_mosaic (Mosaic *mosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info, FlatCorrectionTable *flatcorr); void *setMmos_worker (void *data); -int setMmos_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr); - -void SetMmosInfoInit (SetMmosInfo *info, off_t Nmax, int allocLists, int PoorImages) { +int setMmos_threaded (Catalog *catalog, FlatCorrectionTable *flatcorr); + +void SetMmosInfoInit (SetMmosInfo *info, off_t Nmax, int allocLists) { info->Nfew = 0; info->Nbad = 0; @@ -1014,5 +1014,4 @@ info->Nmax = Nmax; - info->PoorImages = PoorImages; if (allocLists) { @@ -1069,5 +1068,5 @@ } -int setMmos (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { +int setMmos (Catalog *catalog, FlatCorrectionTable *flatcorr) { off_t i, N, Nmax; @@ -1082,5 +1081,5 @@ // so do not run setMmos in threaded mode if PLOTSTUFF is set if (NTHREADS && !PLOTSTUFF) { - int status = setMmos_threaded (catalog, PoorImages, flatcorr); + int status = setMmos_threaded (catalog, flatcorr); return status; } @@ -1089,10 +1088,4 @@ fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE); - - if (PoorImages) { - // XXX use bad stars and measurements for PoorImages? or not? - // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0; - IMAGE_BAD = 0; - } Nmax = 0; @@ -1102,5 +1095,5 @@ SetMmosInfo info; - SetMmosInfoInit (&info, Nmax, TRUE, PoorImages); + SetMmosInfoInit (&info, Nmax, TRUE); // int savePlotDelay = PLOTDELAY; @@ -1121,9 +1114,4 @@ fprintf (stderr, "%d mosaics marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", info.Nfew, info.Nbad, info.Ncal, info.Ngrid, info.Nrel, info.Nsys); - if (PoorImages) { - IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP; - STAR_BAD = ID_OBJ_POOR | ID_OBJ_FEW; - MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA; - } return (TRUE); } @@ -1134,10 +1122,4 @@ off_t j, NimageReal; - /* on PoorImages run, skip good images */ - if (info->PoorImages) { - int bad = myMosaic[0].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP); - if (!bad) return TRUE; - } - FitDataSet *psfStars = &info->psfStars; FitDataSet *kronStars = &info->kronStars; @@ -1147,29 +1129,38 @@ assert (Nmos < Nmosaic); + int Nsecfilt = GetPhotcodeNsecfilt (); + + int minUbercalDist = 1000; + + // unset this flag at start (set below if zeropoint is calculated) + myMosaic->flags &= ~ID_IMAGE_MOSAIC_PHOTCAL; // unset this flag + + // calculate the statistics for both good and bad exposures, but only set Mmos if the rules allow + // we only skip the statistics for nights with too few measurements or exposures + + // if we are fitting TGroup zero points, when we identify the bad TGroups, we proceed + // to fit the mosaics which are from the bad TGroups + + // if we are fitting Mosaic zero points, when we identify the bad nights, we proceed + // to fit the images which are from the bad Mosaics + int badNight = (myMosaic[0].flags & ID_IMAGE_NIGHT_POOR); int badMosaic = (myMosaic[0].flags & ID_IMAGE_MOSAIC_POOR); - - myMosaic->flags &= ~ID_IMAGE_MOSAIC_PHOTCAL; // unset this flag - - // if we are fitting TGroup zero points, when we identify the bad TGroups, we proceed - // to fit the mosaics which are from the bad TGroups - - // if we are fitting Mosaic zero points, when we identify the bad nights, we proceed - // to fit the images which are from the bad Mosaics - - // in BAD_NIGHT mode, we fit ONLY mosaics in bad nights + int useMmos = TRUE; + + // in BAD_NIGHT mode, we fit ONLY mosaics in bad nights (do not fit mosaics from good nights) if (MOSAIC_ZPT_MODE == MOSAIC_ZPT_MODE_BAD_NIGHT) { - if (!badNight) return TRUE; + if (!badNight) useMmos = FALSE; } // in BAD_NIGHT_GOOD_MOSAIC mode, we fit ONLY good mosaics in bad nights if ((MOSAIC_ZPT_MODE == MOSAIC_ZPT_MODE_BAD_NIGHT_GOOD_MOSAIC)) { - if (!badNight) return TRUE; - if (badMosaic) return TRUE; + if (!badNight) useMmos = FALSE; // do not fit good nights + if (badMosaic) useMmos = FALSE; // do not fit bad mosaics } // in GOOD_MOSAIC mode, we fit good mosaics ignoring night state if ((MOSAIC_ZPT_MODE == MOSAIC_ZPT_MODE_GOOD_MOSAIC)) { - if (badMosaic) return TRUE; + if (badMosaic) useMmos = FALSE; // do not fit bad mosaics } @@ -1189,10 +1180,7 @@ // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist); } + // XXX if KEEP_UBERCAL, calculate stats but do not change Mmos if (KEEP_UBERCAL) return TRUE; } - - int Nsecfilt = GetPhotcodeNsecfilt (); - - int minUbercalDist = 1000; // do not modify the calibration for mosaics with partial images loaded (skipCal TRUE) @@ -1204,15 +1192,4 @@ int testImage = FALSE; // testImage |= (abs(myMosaic[0].start - 1324104046) < 10); - // testImage |= (abs(myMosaic[0].start - 1324103823) < 10); - // testImage |= (abs(myMosaic[0].start - 1323003245) < 10); - // testImage |= (abs(myMosaic[0].start - 1323003069) < 10); - // testImage |= (abs(myMosaic[0].start - 1323003125) < 10); - // testImage |= (abs(myMosaic[0].start - 1323003300) < 10); - // testImage |= (abs(myMosaic[0].start - 1323003365) < 10); - // testImage |= (abs(myMosaic[0].start - 1323003191) < 10); - // testImage |= (abs(myMosaic[0].start - 1323003014) < 10); - // testImage |= (abs(myMosaic[0].start - 1323003484) < 10); - // testImage |= (abs(myMosaic[0].start - 1323003419) < 10); - // testImage |= (abs(myMosaic[0].start - 1323002949) < 10); FILE *fout = NULL; @@ -1221,4 +1198,6 @@ snprintf (filename, 64, "test.%05d.%02d.dat", (int) Nmos, npass_output); fout = fopen (filename, "w"); + save_test_mosaic_measures (fout, Nmos, catalog, flatcorr); + fclose (fout); } @@ -1232,22 +1211,4 @@ off_t c = MosaicToCatalog[Nmos][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); - // float Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass, NULL); - 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; - 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, MsysPSF, MrelPSF, Mcal, Mgrid, Mflat, delta, isBad); - } - if (catalog[c].measureT[m].dbFlags & MEAS_BAD) { info->Nbad ++; @@ -1329,38 +1290,63 @@ /* N_onMosaic[Nmos] is all measurements, N is good measurements */ - if (fout) { fclose (fout); } - - /* too few good measurements or too many bad measurements (skip in PoorImages run) */ - if (!info->PoorImages) { - int mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*N_onMosaic[Nmos]); - if (mark) { - if (VERBOSE2) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), (%d < %d) || (%d < %f*"OFF_T_FMT")\n", image[MosaicToImage[Nmos][0]].name, Nmos, N, IMAGE_TOOFEW, N, IMAGE_GOOD_FRACTION, N_onMosaic[Nmos]); } - myMosaic[0].flags |= ID_IMAGE_PHOTOM_FEW; - info->Nfew ++; - if (testImage) { - fprintf (stderr, "NOTE: *** marked test image poor : %d %d %d***\n", (int) N, (int) IMAGE_TOOFEW, (int) (IMAGE_GOOD_FRACTION*N_onMosaic[Nmos])); - } - return TRUE; // skip mosaics with too few good measurements - } else { - myMosaic[0].flags &= ~ID_IMAGE_PHOTOM_FEW; - } + /* skip mosaics with too few good measurements or too many bad measurements */ + int mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*N_onMosaic[Nmos]); + if (mark) { + myMosaic->flags |= ID_IMAGE_PHOTOM_FEW; + myMosaic->McalPSF = 0.0; + myMosaic->McalAPER = 0.0; + myMosaic->dMcal = NAN; + myMosaic->stdev = NAN; + myMosaic->dMmin = NAN; + myMosaic->dMmax = NAN; + myMosaic->McalChiSq = NAN; + myMosaic->nFitPhotom = 0; + myMosaic->nValPhotom = N; + info->Nfew ++; + if (testImage) { + fprintf (stderr, "NOTE: *** marked test image poor : %d %d %d***\n", (int) N, (int) IMAGE_TOOFEW, (int) (IMAGE_GOOD_FRACTION*N_onMosaic[Nmos])); + } + return TRUE; // skip mosaics with too few good measurements + } else { + myMosaic[0].flags &= ~ID_IMAGE_PHOTOM_FEW; + } + + 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; } fit1d_irls (psfStars, N); fit1d_irls (brightStars, Nbright); + fit1d_irls (kronStars, N); + + if (useMmos) { + myMosaic->McalPSF = psfStars->bSaveArray[0][0]; + myMosaic->McalAPER = kronStars->bSaveArray[0][0]; + myMosaic->flags |= ID_IMAGE_MOSAIC_PHOTCAL; // set this flag + } else { + myMosaic->McalPSF = 0.0; + myMosaic->McalAPER = 0.0; + } // for now, I have no reason to measure these separately for camera-level images - 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; - myMosaic[0].nValPhotom = N; - - fit1d_irls (kronStars, N); - myMosaic->McalAPER = kronStars->bSaveArray[0][0]; + myMosaic->dMcal = psfStars->bSigma[0]; + myMosaic->stdev = psfStars->sigma; + myMosaic->McalChiSq = psfStars->chisq; + myMosaic->nFitPhotom = psfStars->Nmeas; + myMosaic->nValPhotom = N; + + // bright end scatter + myMosaic->dMsys = brightStars->sigma; // stdev of bright star mags if (testImage) { fprintf (stderr, "test image %d (%d) %f %f %d ... ", (int) Nmos, myMosaic->start, myMosaic->McalPSF, myMosaic->dMcal, myMosaic->nFitPhotom); + fprintf (stderr, "%f %f : %f\n", myMosaic[0].McalPSF, myMosaic[0].dMsys, myMosaic[0].McalChiSq); } @@ -1368,18 +1354,4 @@ fprintf (stderr, "Mmos: %6.3f +/- %6.3f %5d %5d | %s\n", myMosaic->McalPSF, myMosaic->dMcal, myMosaic->nFitPhotom, N, image[MosaicToImage[Nmos][0]].name); plot_setMcal (psfStars->alldata->yVector, N); - } - - // bright end scatter - fit1d_irls (brightStars, Nbright); - myMosaic->dMsys = brightStars->sigma; // stdev of bright star mags - - if (myMosaic[0].McalPSF < -CLOUD_TOLERANCE) { - myMosaic[0].McalPSF = 0.0; - } - - myMosaic->flags |= ID_IMAGE_MOSAIC_PHOTCAL; // set this flag - - if (testImage) { - fprintf (stderr, "%f %f : %f\n", myMosaic[0].McalPSF, myMosaic[0].dMsys, myMosaic[0].McalChiSq); } @@ -1400,5 +1372,33 @@ } -int setMmos_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { +int save_test_mosaic_measures (FILE *fout, int Nmos, Catalog *catalog, FlatCorrectionTable *flatcorr) { + + int Nsecfilt = GetPhotcodeNsecfilt (); + + for (int j = 0; j < N_onMosaic[Nmos]; j++) { + + 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 + float Mcal = getMcal (m, c, MAG_CLASS_PSF); + // float Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass, NULL); + 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; + 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, MsysPSF, MrelPSF, Mcal, Mgrid, Mflat, delta, isBad); + } + return TRUE; +} + +int setMmos_threaded (Catalog *catalog, FlatCorrectionTable *flatcorr) { int i; @@ -1408,10 +1408,4 @@ fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE); - - if (PoorImages) { - // XXX use bad stars and measurements for PoorImages? or not? - // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0; - IMAGE_BAD = 0; - } off_t Nmax = 0; @@ -1421,5 +1415,5 @@ SetMmosInfo summary; - SetMmosInfoInit (&summary, Nmax, FALSE, PoorImages); + SetMmosInfoInit (&summary, Nmax, FALSE); pthread_attr_t attr; @@ -1445,7 +1439,7 @@ threadinfo[i].flatcorr = flatcorr; - // we do NOT allocate the arrays here, we only supply basic info (Nmax, Nloop, - // PoorImages) used in the threads to allocate the arrays and set the MaxIterations - SetMmosInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages); + // we do NOT allocate the arrays here, we only supply basic info (Nmax, Nloop) used in + // the threads to allocate the arrays and set the MaxIterations + SetMmosInfoInit (&threadinfo[i].info, Nmax, FALSE); pthread_create (&threads[i], NULL, setMmos_worker, &threadinfo[i]); } @@ -1492,9 +1486,4 @@ npass_output ++; - if (PoorImages) { - IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP; - STAR_BAD = ID_OBJ_POOR | ID_OBJ_FEW; - MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA; - } return TRUE; } @@ -1505,5 +1494,5 @@ SetMmosInfo results; - SetMmosInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages); + SetMmosInfoInit (&results, threadinfo->info.Nmax, TRUE); while (1) { @@ -1892,5 +1881,4 @@ 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; @@ -1906,11 +1894,6 @@ 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); @@ -1927,5 +1910,5 @@ int mark = FALSE; - if (mosaic[i].flags & (ID_IMAGE_PHOTOM_FEW)) { + if (mosaic[i].flags & ID_IMAGE_PHOTOM_FEW) { mark = TRUE; Nfew ++; Index: /trunk/Ohana/src/relphot/src/StarOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/StarOps.c (revision 41555) +++ /trunk/Ohana/src/relphot/src/StarOps.c (revision 41556) @@ -662,4 +662,6 @@ // this function only operations on the PSF magnitudes + +// XXX EAM clean_measures is no longer used -- we rely on IRLS to down-weight outliers # define NSIGMA_CLIP 3.0 Index: /trunk/Ohana/src/relphot/src/TGroupOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/TGroupOps.c (revision 41555) +++ /trunk/Ohana/src/relphot/src/TGroupOps.c (revision 41556) @@ -194,10 +194,4 @@ for (off_t i = 0; i < Nsubset; i++) { ImageToTGroup[i] = NULL; - - // ignore non-GPC1 images - // if (!isGPC1chip(subset[i].photcode)) { - // Nsimple ++; - // continue; - // } // ignore mosaic images (photcode == 0) @@ -459,4 +453,7 @@ // this tgroup has been identified as poor. tgroup->flags |= ID_IMAGE_NIGHT_POOR; + tgroup->McalPSF = 0.0; + tgroup->McalAPER = 0.0; + tgroup->dKlam = 0.0; // all associated mosaics should be marked as coming from a bad night @@ -473,8 +470,4 @@ image[im].flags |= ID_IMAGE_NIGHT_POOR; } - - tgroup->McalPSF = 0.0; - tgroup->McalAPER = 0.0; - tgroup->dKlam = 0.0; } @@ -552,11 +545,4 @@ fprintf (stderr, "Matched %d detections to tgroups\n", Nmatch); - // XXX print results for testing - for (int i = 0; FALSE && (i < NtgroupTimes); i++) { - TGroup *tgroup = tgroupTimes[i][0].byCode; - for (off_t j = 0; j < tgroupTimes[i][0].nCode; j++) { - if (0) fprintf (stderr, " TGROUP time %.0f photcode %d, Nmeasure: %d, Nimage: %d\n", ohana_sec_to_mjd(tgroupTimes[i][0].start), tgroup[j].photcode, (int) tgroup[j].Nmeasure, (int) tgroup[j].Nimage); - } - } return (TRUE); } @@ -574,24 +560,5 @@ TGroup *myGroup = ImageToTGroup[idx]; - TGTimes *parent = (TGTimes *) myGroup->parent; - if (!myGroup) { - // Image *image = getimage(idx); - // fprintf (stderr, "unmatched image %s\n", image[0].name); - // skip measurements from simple chips (not matched to a tgroup by definition) - return; - } - - // test to check we got the right match: - { - Image *image = getimage(idx); - // XXX we are now matching with just tzero. be careful for cameras with drift - // unsigned int imageStart = image[0].tzero - MAX(0.01*image[0].trate*image[0].NY, 1); - unsigned int imageStart = image[0].tzero; - // XXX to do this check I need a pointer to the group master - if ((imageStart < parent->start) || (imageStart > parent->stop)) { - fprintf (stderr, "error in image to tgroup match\n"); - abort(); - } - } + if (!myGroup) return; // this measurement is on one of my tgroups, mark it as mine. @@ -611,5 +578,5 @@ } -// XXX return the error as well as the value +// dZpt (if supplied) returns the error on Mgrp float getMgrp (off_t meas, int cat, float airmass, float *dZpt) { @@ -638,6 +605,4 @@ } -// XXX remove the lists and replace with StatDataSet structs? -// XXX move the StatDataSetAlloc up? typedef struct { int Nfew; @@ -650,5 +615,4 @@ int Nskip; off_t Nmax; - int PoorImages; FitDataSet psfStars; FitDataSet kronStars; @@ -669,7 +633,7 @@ int setMgrp_tgroup (TGroup *tgroup, off_t Nmos, Image *image, Catalog *catalog, SetMgrpInfo *info, FlatCorrectionTable *flatcorr); void *setMgrp_worker (void *data); -int setMgrp_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr); - -void SetMgrpInfoInit (SetMgrpInfo *info, off_t Nmax, int allocLists, int PoorImages) { +int setMgrp_threaded (Catalog *catalog, FlatCorrectionTable *flatcorr); + +void SetMgrpInfoInit (SetMgrpInfo *info, off_t Nmax, int allocLists) { info->Nfew = 0; info->Nbad = 0; @@ -682,5 +646,4 @@ info->Nmax = Nmax; - info->PoorImages = PoorImages; if (allocLists) { @@ -744,5 +707,5 @@ } -int setMgrp (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { +int setMgrp (Catalog *catalog, FlatCorrectionTable *flatcorr) { off_t N; @@ -755,5 +718,5 @@ // so do not run setMgrp in threaded mode if PLOTSTUFF is set if (NTHREADS && !PLOTSTUFF) { - int status = setMgrp_threaded (catalog, PoorImages, flatcorr); + int status = setMgrp_threaded (catalog, flatcorr); return status; } @@ -770,5 +733,5 @@ SetMgrpInfo info; - SetMgrpInfoInit (&info, Nmax, TRUE, PoorImages); + SetMgrpInfoInit (&info, Nmax, TRUE); for (off_t i = 0; i < NtgroupTimes; i++) { @@ -776,5 +739,6 @@ for (off_t j = 0; j < tgroupTimes[i][0].nCode; j++) { setMgrp_tgroup (&tgroup[j], i, image, catalog, &info, flatcorr); - // fprintf (stderr, "TGROUP time %.0f photcode %d Mcal: %f, dK: %f, dMcal: %f, sigma: %f, chisq: %f, %d of %d, limiting negative clouds to %f\n", ohana_sec_to_mjd(tgroupTimes[i][0].start), tgroup[j].photcode, tgroup[j].McalPSF, tgroup[j].dKlam, tgroup[j].dMcal, tgroup[j].stdev, tgroup[j].McalChiSq, tgroup[j].nFitPhotom, tgroup[j].nValPhotom, CLOUD_TOLERANCE); + // fprintf (stderr, "TGROUP time %.0f photcode %d Mcal: %f, dK: %f, dMcal: %f, sigma: %f, chisq: %f, %d of %d, limiting negative clouds to %f\n", + // ohana_sec_to_mjd(tgroupTimes[i][0].start), tgroup[j].photcode, tgroup[j].McalPSF, tgroup[j].dKlam, tgroup[j].dMcal, tgroup[j].stdev, tgroup[j].McalChiSq, tgroup[j].nFitPhotom, tgroup[j].nValPhotom, CLOUD_TOLERANCE); } } @@ -799,12 +763,8 @@ int setMgrp_tgroup (TGroup *myTGroup, off_t Ngrp, Image *image, Catalog *catalog, SetMgrpInfo *info, FlatCorrectionTable *flatcorr) { - off_t j, NimageReal; - - // At some point, we identify good and bad tgroups. We mark bad tgroups and do not fit - // a zero point for them. Instead we will fit those exposures or images separately - if (TGROUP_ZPT_MODE == TGROUP_ZPT_MODE_GOOD_NIGHT) { - int bad = myTGroup->flags & ID_IMAGE_NIGHT_POOR; - if (bad) return TRUE; - } + off_t j; + + // calculate the statistics for both good and bad nights, but only set Mcal for the good nights + // we only skip the statistics for nights with too few measurements or exposures FitDataSet *psfStars = &info->psfStars; @@ -815,10 +775,6 @@ assert (Ngrp < NtgroupTimes); - // Image *imageReal = getimages (&NimageReal, NULL); returned pointer is not used - getimages (&NimageReal, NULL); - int Nsecfilt = GetPhotcodeNsecfilt (); - int testImage = FALSE; TGTimes *tgroup = (TGTimes *) myTGroup->parent; double mjdStart = ohana_sec_to_mjd (tgroup->start); @@ -827,14 +783,15 @@ 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); + // in clean_tgroups (run on each iteration), we identify good and bad tgroups. We mark + // bad tgroups and do not save a zero point for them. Instead we will fit those + // exposures or images separately + int useMgrp = TRUE; + if (TGROUP_ZPT_MODE == TGROUP_ZPT_MODE_GOOD_NIGHT) { + if (myTGroup->flags & ID_IMAGE_NIGHT_POOR) useMgrp = FALSE; + } + + // for testing, supply the MJD of a night or nights here to dump the full list of measurements + int testImage = FALSE; + // testImage |= (abs(mjdStart - 57353) < 0.1); FILE *fout = NULL; @@ -905,5 +862,5 @@ float MsysKron = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_KRON); - // XXX byCode : here we are only checking for equiv photcodes (not active photcodes) + // here we are only checking for equiv photcodes (not active photcodes) PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode); if (!code) goto skip; @@ -942,15 +899,18 @@ } - /* too few good measurements or too many bad measurements (skip in PoorImages run) */ + /* skip nights with too few good measurements or too many bad measurements */ int mark = (N < NIGHT_TOOFEW) || (N < NIGHT_GOOD_FRACTION*myTGroup->Nmeasure); if (mark) { - if (TRUE) { - TGTimes *mygroup = (TGTimes *) myTGroup->parent; - 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->McalPSF = 0.0; + myTGroup->McalAPER = 0.0; + myTGroup->dMcal = NAN; + myTGroup->stdev = NAN; + myTGroup->dMmin = NAN; + myTGroup->dMmax = NAN; + myTGroup->McalChiSq = NAN; + myTGroup->nFitPhotom = 0; myTGroup->nValPhotom = N; + myTGroup->dKlam = (TGROUP_FIT_AIRMASS) ? psfStars->bSaveArray[1][0] : 0; info->Nfew ++; if (testImage) { @@ -986,9 +946,17 @@ fit1d_irls (psfStars, N); - - // if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mgrp: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N); - - // for now, I have no reason to measure these separately for camera-level images - myTGroup->McalPSF = psfStars->bSaveArray[0][0]; + fit1d_irls (kronStars, N); + fit1d_irls (brightStars, Nbright); + + if (useMgrp) { + myTGroup->McalPSF = psfStars->bSaveArray[0][0]; + myTGroup->McalAPER = kronStars->bSaveArray[0][0]; // XXX this does not make sense: I need to apply the airmass slope calculated above first + myTGroup->flags |= ID_IMAGE_TGROUP_PHOTCAL; // set this flag (unset by default) + } else { + myTGroup->McalPSF = 0.0; + myTGroup->McalAPER = 0.0; + } + + // record statistics on the fit regardless if it is used or not myTGroup->dMcal = psfStars->bSigma[0]; myTGroup->stdev = psfStars->sigma; @@ -1000,32 +968,15 @@ 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 - fit1d_irls (kronStars, N); - myTGroup->McalAPER = kronStars->bSaveArray[0][0]; + // bright end scatter (systematic error) + myTGroup->dMsys = brightStars->sigma; if (testImage) { TGTimes *parent = (TGTimes *) myTGroup->parent; fprintf (stderr, "test night %.0f aper: %f %f %d ... ", ohana_sec_to_mjd(parent->start), myTGroup->McalPSF, myTGroup->dMcal, myTGroup->nFitPhotom); + fprintf (stderr, "%f %f : %f\n", myTGroup->McalPSF, myTGroup->dMsys, myTGroup->McalChiSq); } if (PLOTSTUFF) { fprintf (stderr, "Mgrp: %6.3f +/- %6.3f %5d of %5d | %.0f\n", myTGroup->McalPSF, myTGroup->dMcal, myTGroup->nFitPhotom, N, ohana_sec_to_mjd(tgroup->start)); plot_setMcal (psfStars->alldata->yVector, N); - } - - // bright end scatter - // XXX this does not make sense: I need to apply the airmass slope calculated above first - // redo this by calculating the corrected bright stars mags - fit1d_irls (brightStars, Nbright); - myTGroup->dMsys = brightStars->sigma; - - // keep this?? - if (myTGroup->McalPSF < -CLOUD_TOLERANCE) { - myTGroup->McalPSF = 0.0; - } - - myTGroup->flags |= ID_IMAGE_TGROUP_PHOTCAL; // set this flag - - if (testImage) { - fprintf (stderr, "%f %f : %f\n", myTGroup->McalPSF, myTGroup->dMsys, myTGroup->McalChiSq); } @@ -1070,5 +1021,5 @@ } -int setMgrp_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { +int setMgrp_threaded (Catalog *catalog, FlatCorrectionTable *flatcorr) { int i; @@ -1077,6 +1028,4 @@ Image *image = getimages (&N, NULL); - fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE); - off_t Nmax = 0; for (off_t i = 0; i < NtgroupTimes; i++) { @@ -1088,5 +1037,5 @@ SetMgrpInfo summary; - SetMgrpInfoInit (&summary, Nmax, FALSE, PoorImages); + SetMgrpInfoInit (&summary, Nmax, FALSE); pthread_attr_t attr; @@ -1112,7 +1061,7 @@ threadinfo[i].flatcorr = flatcorr; - // we do NOT allocate the arrays here, we only supply basic info (Nmax, - // PoorImages) used in the threads to allocate the arrays and set the MaxIterations - SetMgrpInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages); + // we do NOT allocate the arrays here, we only supply basic info used in the threads + // to allocate the arrays and set the MaxIterations + SetMgrpInfoInit (&threadinfo[i].info, Nmax, FALSE); pthread_create (&threads[i], NULL, setMgrp_worker, &threadinfo[i]); } @@ -1169,5 +1118,5 @@ SetMgrpInfo results; - SetMgrpInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages); // allocate list, dlist arrays here + SetMgrpInfoInit (&results, threadinfo->info.Nmax, TRUE); // allocate list, dlist arrays here while (1) { @@ -1194,6 +1143,4 @@ } -// XXX should I split these stats out by photcode? -// XXX double check the IMAGE_BAD filter: StatType statsTGroupM (Catalog *catalog) { OHANA_UNUSED_PARAM(catalog); @@ -1209,6 +1156,5 @@ ALLOCATE (dlist, double, NtgroupTimes*Nphotcodes); - // is ID_IMAGE_PHOTOM_POOR used for tgroups? - int NIGHT_BAD = ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR; + int NIGHT_BAD = ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW; int n = 0; @@ -1316,12 +1262,13 @@ ALLOCATE (slist, double, NtgroupTimes*Nphotcodes); + // measure the good/bad statistics using nights which are actually calibrated + // but they will be applied to all nights below (allowing bad nights to become good) int N = 0; for (int i = 0; i < NtgroupTimes; i++) { TGroup *tgroup = tgroupTimes[i][0].byCode; for (int j = 0; j < tgroupTimes[i][0].nCode; j++) { - if (tgroup[j].flags & IMAGE_BAD) continue; + if (!(tgroup[j].flags & ID_IMAGE_TGROUP_PHOTCAL)) continue; mlist[N] = tgroup[j].McalChiSq; slist[N] = tgroup[j].stdev; // stdev of all measurements - // slist[N] = tgroup[j].dMsys; // stdev of bright measurements N++; } @@ -1491,4 +1438,5 @@ } +// temporary test output void dump_tgroups (Catalog *catalog, int Npass) { @@ -1606,2 +1554,9 @@ fclose (fout); } + +/* + for testing, put lines like these in relphot_images.c within the iteration loop + dump_tgroups (catalog, i + 40); + dump_catalog (catalog, 0, i + 40); // for a test, just dump a specific catalog + dump_tgroup_imstats (i + 40); +*/ Index: /trunk/Ohana/src/relphot/src/relphot_images.c =================================================================== --- /trunk/Ohana/src/relphot/src/relphot_images.c (revision 41555) +++ /trunk/Ohana/src/relphot/src/relphot_images.c (revision 41556) @@ -41,6 +41,4 @@ load_images (&db, skylist, &UserPatch, FALSE, USE_ALL_IMAGES); MARKTIME("-- load images: %f sec\n", dtime); - - // checkImages("load images"); /* unlock, if we can (else, unlocked below) */ @@ -63,6 +61,4 @@ /* add in a loop over the catalogs calling dvo_catalog_chipcoords */ - // checkImages("load catalog"); - /* match measurements with images, mosaics */ initImageBins (catalog, Ncatalog, TRUE); @@ -87,5 +83,4 @@ setExclusions (catalog, Ncatalog, TRUE); - // checkImages("set Exclu"); global_stats (catalog, Ncatalog, flatcorr, 0); @@ -118,21 +113,9 @@ 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); - setMmos (catalog, FALSE, flatcorr); - setMgrp (catalog, FALSE, flatcorr); - dump_tgroups (catalog, i + 40); - dump_catalog (catalog, 0, i + 40); // for a test, just dump a specific catalog - dump_tgroup_imstats (i + 40); + + setMcal (catalog, flatcorr); + setMmos (catalog, flatcorr); + setMgrp (catalog, flatcorr); setMgrid (catalog, flatcorr); @@ -166,10 +149,4 @@ if (USE_GRID) dump_grid (); - /* set Mcal & Mmos for bad images */ - int onlyPoorImages = !CALIBRATE_STACKS_AND_WARPS; - onlyPoorImages = FALSE; // XXX deactivate for now - setMcal (catalog, onlyPoorImages, flatcorr); // max loop since this is a final pass - setMmos (catalog, onlyPoorImages, flatcorr); // max loop since this is a final pass - setMgrp (catalog, onlyPoorImages, flatcorr); // max loop since this is a final pass MARKTIME("-- finalize Mcal values: %f sec\n", dtime); @@ -177,67 +154,8 @@ setMcalFromTGroups (); // copy per-tgroup calibrations to the images - if (SAVE_IMAGE_UPDATES) { - - FITS_DB dbX; - off_t *LineNumber; - Image *image, *subset; - off_t Nimage, Nsubset, i, Nx; - char filename[1024]; - - gfits_db_init (&dbX); - dbX.lockstate = LCK_XCLD; - dbX.timeout = 60.0; - dbX.mode = db.mode; - dbX.format = db.format; - - snprintf (filename, 1024, "%s.bck", ImageCat); - if (!gfits_db_lock (&dbX, filename)) { - fprintf (stderr, "can't lock backup image image catalog\n"); - return (FALSE); - } - - // apply the changes from the image subset to the full table: - - // convert database table to internal structure (binary to Image) - // 'image' points to the same memory as db->ftable->buffer - image = gfits_table_get_Image (&db.ftable, &Nimage, &db.scaledValue, &db.nativeOrder); - if (!image) { - fprintf (stderr, "ERROR: failed to read images\n"); - exit (2); - } - gfits_scan (db.ftable.header, "NAXIS1", OFF_T_FMT, 1, &Nx); - - subset = getimages (&Nsubset, &LineNumber); - for (i = 0; i < Nsubset; i++) { - if (LineNumber[i] == -1) continue; - memcpy (&image[LineNumber[i]], &subset[i], Nx); - } - - // copy Images.dat data (all or only the subset / vtable elements?) I think I need to dump - // the entire Image table, but with the updates in place I think this says: re-work the - // ftable/vtable usage in this program to be more sensible... - gfits_copy_header (&db.header, &dbX.header); - gfits_copy_matrix (&db.matrix, &dbX.matrix); - gfits_copy_header (&db.theader, &dbX.theader); - gfits_copy_ftable (&db.ftable, &dbX.ftable); - - dbX.ftable.header = &dbX.theader; - dbX.virtual = FALSE; - - // save Images.dat using the copied structure - if (UPDATE_CATFORMAT) { - // 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); - SavePhotcodesFITS (photcodeFile); - } - dvo_image_save (&dbX, VERBOSE); - dvo_image_unlock (&dbX); - gfits_db_free (&dbX); - MARKTIME("-- save Image.dat.bck: %f sec\n", dtime); - } + // calculate and save per-chip residuals here: + MagResidSave ("mag.resid.fits", catalog); + + save_images_updates (&db); /* at this point, we have correct cal coeffs in the image/mosaic structures */ Index: /trunk/Ohana/src/relphot/src/relphot_parallel_images.c =================================================================== --- /trunk/Ohana/src/relphot/src/relphot_parallel_images.c (revision 41555) +++ /trunk/Ohana/src/relphot/src/relphot_parallel_images.c (revision 41556) @@ -112,6 +112,6 @@ // set the image (Mcal) and mosaic (Mmos) zero point offsets given the mean mags and measurements - setMcal (catalog, FALSE, flatcorr); - setMmos (catalog, FALSE, flatcorr); + setMcal (catalog, flatcorr); + setMmos (catalog, flatcorr); // setMgrp (catalog, FALSE, flatcorr); XXX think this through: this may not make sense for tgroups MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime); @@ -144,9 +144,4 @@ client_logger_message ("all hosts are done the loops\n"); - /* 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); - // setMgrp (catalog, onlyPoorImages, flatcorr); XXX see above MARKTIME("-- finalize Mcal values: %f sec\n", dtime); Index: /trunk/Ohana/src/relphot/src/relphot_parallel_regions.c =================================================================== --- /trunk/Ohana/src/relphot/src/relphot_parallel_regions.c (revision 41555) +++ /trunk/Ohana/src/relphot/src/relphot_parallel_regions.c (revision 41556) @@ -1,4 +1,3 @@ # include "relphot.h" -int save_images_backup (FITS_DB *db); int relphot_parallel_regions (SkyTable *sky) { @@ -101,39 +100,2 @@ return TRUE; } - -int save_images_backup (FITS_DB *db) { - - FITS_DB dbX; - char filename[1024]; - - INITTIME; - - gfits_db_init (&dbX); - dbX.lockstate = LCK_XCLD; - dbX.timeout = 60.0; - dbX.mode = db->mode; - dbX.format = db->format; - - snprintf (filename, 1024, "%s.bck", ImageCat); - if (!gfits_db_lock (&dbX, filename)) { - fprintf (stderr, "can't lock backup image image catalog\n"); - return (FALSE); - } - - // copy Images.dat data (all or only the subset / vtable elements?) I think I need to dump - // the entire Image table, but with the updates in place I think this says: re-work the - // ftable/vtable usage in this program to be more sensible... - gfits_copy_header (&db->header, &dbX.header); - gfits_copy_matrix (&db->matrix, &dbX.matrix); - gfits_copy_header (&db->theader, &dbX.theader); - gfits_copy_ftable (&db->ftable, &dbX.ftable); - - dbX.ftable.header = &dbX.theader; - - dvo_image_save (&dbX, VERBOSE); - dvo_image_unlock (&dbX); - MARKTIME("-- save Image.dat.bck: %f sec\n", dtime); - - gfits_db_free (&dbX); - return TRUE; -} Index: /trunk/Ohana/src/relphot/src/save_images.c =================================================================== --- /trunk/Ohana/src/relphot/src/save_images.c (revision 41556) +++ /trunk/Ohana/src/relphot/src/save_images.c (revision 41556) @@ -0,0 +1,108 @@ +# include "relphot.h" + +// XXX there are two versions of this code: +// save_image_updates () is called by relphot_images +// save_image_backup () is called by relphot_parallel_regions + +// rationalize these to a single function... + +int save_images_updates (FITS_DB *db) { + + if (!SAVE_IMAGE_UPDATES) return TRUE; + + FITS_DB dbX; + gfits_db_init (&dbX); + dbX.lockstate = LCK_XCLD; + dbX.timeout = 60.0; + dbX.mode = db->mode; + dbX.format = db->format; + + char filename[1024]; + snprintf (filename, 1024, "%s.bck", ImageCat); // ImageCat is global + if (!gfits_db_lock (&dbX, filename)) { + fprintf (stderr, "can't lock backup image image catalog\n"); + return FALSE; + } + + // apply the changes from the image subset to the full table: + + // convert database table to internal structure (binary to Image) + // 'image' points to the same memory as db->ftable->buffer + + off_t Nimage, Nx; + Image *image = gfits_table_get_Image (&db->ftable, &Nimage, &db->scaledValue, &db->nativeOrder); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } + gfits_scan (db->ftable.header, "NAXIS1", OFF_T_FMT, 1, &Nx); + + off_t *LineNumber; + off_t Nsubset; + Image *subset = getimages (&Nsubset, &LineNumber); + for (off_t i = 0; i < Nsubset; i++) { + if (LineNumber[i] == -1) continue; + memcpy (&image[LineNumber[i]], &subset[i], Nx); + } + + // copy Images.dat data (all or only the subset / vtable elements?) I think I need to dump + // the entire Image table, but with the updates in place I think this says: re-work the + // ftable/vtable usage in this program to be more sensible... + gfits_copy_header (&db->header, &dbX.header); + gfits_copy_matrix (&db->matrix, &dbX.matrix); + gfits_copy_header (&db->theader, &dbX.theader); + gfits_copy_ftable (&db->ftable, &dbX.ftable); + + dbX.ftable.header = &dbX.theader; + dbX.virtual = FALSE; + + // save Images.dat using the copied structure + if (UPDATE_CATFORMAT) { + // 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); + SavePhotcodesFITS (photcodeFile); + } + + dvo_image_save (&dbX, VERBOSE); + dvo_image_unlock (&dbX); + gfits_db_free (&dbX); + return TRUE; +} + +int save_images_backup (FITS_DB *db) { + + FITS_DB dbX; + char filename[1024]; + + gfits_db_init (&dbX); + dbX.lockstate = LCK_XCLD; + dbX.timeout = 60.0; + dbX.mode = db->mode; + dbX.format = db->format; + + snprintf (filename, 1024, "%s.bck", ImageCat); + if (!gfits_db_lock (&dbX, filename)) { + fprintf (stderr, "can't lock backup image image catalog\n"); + return (FALSE); + } + + // copy Images.dat data (all or only the subset / vtable elements?) I think I need to dump + // the entire Image table, but with the updates in place I think this says: re-work the + // ftable/vtable usage in this program to be more sensible... + gfits_copy_header (&db->header, &dbX.header); + gfits_copy_matrix (&db->matrix, &dbX.matrix); + gfits_copy_header (&db->theader, &dbX.theader); + gfits_copy_ftable (&db->ftable, &dbX.ftable); + + dbX.ftable.header = &dbX.theader; + + dvo_image_save (&dbX, VERBOSE); + dvo_image_unlock (&dbX); + + gfits_db_free (&dbX); + return TRUE; +}