Index: trunk/Ohana/src/relphot/Makefile =================================================================== --- trunk/Ohana/src/relphot/Makefile (revision 41461) +++ trunk/Ohana/src/relphot/Makefile (revision 41462) @@ -77,4 +77,5 @@ $(SRC)/indexCatalog.$(ARCH).o \ $(SRC)/client_logger.$(ARCH).o \ +$(SRC)/ZeroPointModes.$(ARCH).o \ $(SRC)/setExclusions.$(ARCH).o \ $(SRC)/setMrelFinal.$(ARCH).o \ @@ -124,4 +125,5 @@ $(SRC)/synthetic_zpts.$(ARCH).o \ $(SRC)/client_logger.$(ARCH).o \ +$(SRC)/ZeroPointModes.$(ARCH).o \ $(SRC)/setExclusions.$(ARCH).o \ $(SRC)/setMrelFinal.$(ARCH).o \ @@ -143,4 +145,5 @@ TEST_FITIRLS = \ $(SRC)/fit1d_irls.$(ARCH).o \ +$(SRC)/extra.$(ARCH).o \ $(SRC)/test_fitirls.$(ARCH).o Index: trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- trunk/Ohana/src/relphot/include/relphot.h (revision 41461) +++ trunk/Ohana/src/relphot/include/relphot.h (revision 41462) @@ -84,7 +84,9 @@ float secz; float ubercalDist; - unsigned short nFitPhotom; + unsigned int nFitPhotom; + unsigned int nValPhotom; unsigned int flags; char skipCal; // if TRUE, this mosaic is incomplete and should not be calibrated + char inTGroup; Coords coords; } Mosaic; @@ -97,6 +99,8 @@ float dMcal; float dMsys; + float stdev; float McalChiSq; - unsigned short nFitPhotom; + unsigned int nFitPhotom; + unsigned int nValPhotom; unsigned int flags; @@ -109,4 +113,8 @@ off_t *measure; off_t *catalog; + + off_t Nmosaic; + off_t NMOSAIC; + Mosaic **mosaic; // pointer to the mosaic structures void *parent; @@ -185,4 +193,5 @@ int Nlist; // total number of measurements in list (may be more than the number to use) int Nbootstrap; // number of bootstrap iterations + int MaxIterations; // number of IRLS iterations // input data & parameters: @@ -387,6 +396,11 @@ int VERBOSE2; int VERBOSE_IMAGE; + int MOSAIC_ZEROPT; int TGROUP_ZEROPT; + +int FREEZE_IMAGES; +int FREEZE_MOSAICS; + int TGROUP_FIT_AIRMASS; @@ -414,6 +428,4 @@ double IMAGE_GOOD_FRACTION; int IMAGE_BAD; -int FREEZE_IMAGES; -int FREEZE_MOSAICS; int CALIBRATE_STACKS_AND_WARPS; int USE_GRID; @@ -461,4 +473,13 @@ SkyRegion UserPatch; char *UserCatalog; + +enum {TGROUP_ZPT_MODE_NONE, TGROUP_ZPT_MODE_GOOD_NIGHT, TGROUP_ZPT_MODE_ALL, }; +int TGROUP_ZPT_MODE; + +enum {MOSAIC_ZPT_MODE_NONE, MOSAIC_ZPT_MODE_BAD_NIGHT, MOSAIC_ZPT_MODE_GOOD_MOSAIC, MOSAIC_ZPT_MODE_BAD_NIGHT_GOOD_MOSAIC, MOSAIC_ZPT_MODE_ALL, }; +int MOSAIC_ZPT_MODE; + +enum {IMAGE_ZPT_MODE_NONE, IMAGE_ZPT_MODE_BAD_NIGHT, IMAGE_ZPT_MODE_BAD_MOSAIC, IMAGE_ZPT_MODE_BAD_NIGHT_BAD_MOSAIC, IMAGE_ZPT_MODE_ALL, }; +int IMAGE_ZPT_MODE; int USE_ALL_IMAGES; @@ -526,5 +547,5 @@ float getMgrid PROTO((off_t meas, int cat)); float getMmos PROTO((off_t meas, int cat)); -float getMgrp PROTO((off_t meas, int cat, float airmass)); +float getMgrp PROTO((off_t meas, int cat, float airmass, float *dZpt)); float getMrel PROTO((Catalog *catalog, off_t meas, int cat, dvoMagClassType class, dvoMagSourceType source)); short getUbercalDist PROTO((off_t meas, int cat)); @@ -601,14 +622,17 @@ int reload_images PROTO((FITS_DB *db)); int setExclusions PROTO((Catalog *catalog, int Ncatalog, int verbose)); -void setMcal PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr)); +void setMgrid PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); void setMcalFromMosaics PROTO((void)); void setMcalFromTGroups PROTO((void)); int setMcalOutput PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); -void setMgrid PROTO((Catalog *catalog, 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, int Nloop, int Poor, FlatCorrectionTable *flatcorr)); +int setMmos PROTO((Catalog *catalog, int Nloop, int Poor, FlatCorrectionTable *flatcorr)); +int setMgrp PROTO((Catalog *catalog, int Nloop, int Poor, FlatCorrectionTable *flatcorr)); + int setMrel PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); void setMrelFinal PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage)); int setMrelOutput PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); + int setMave PROTO((Catalog *catalog, int Ncatalog)); void set_ZP PROTO((double ZERO)); @@ -627,4 +651,8 @@ StatType statsMosaicX PROTO((Catalog *catalog)); StatType statsMosaicdM PROTO((Catalog *catalog)); +StatType statsTGroupM PROTO((Catalog *catalog)); +StatType statsTGroupX PROTO((Catalog *catalog)); +StatType statsTGroupdM PROTO((Catalog *catalog)); + StatType statsStarN PROTO((Catalog *catalog, int Ncatalog, int Nsec, int seccode, FlatCorrectionTable *flatcorr)); StatType statsStarS PROTO((Catalog *catalog, int Ncatalog, int Nsec)); @@ -795,2 +823,7 @@ int save_images_backup (FITS_DB *db); + +void dump_tgroups (Catalog *catalog, int Npass); +void dump_catalog (Catalog *catalog, off_t c, int Npass); + +void SetZeroPointModes (int nloop); Index: trunk/Ohana/src/relphot/src/GridOps.c =================================================================== --- trunk/Ohana/src/relphot/src/GridOps.c (revision 41461) +++ trunk/Ohana/src/relphot/src/GridOps.c (revision 41462) @@ -415,5 +415,5 @@ // skip mosaics marked as BAD - Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass); + Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass, NULL); if (isnan(Mgrp)) { Ngrp ++; @@ -575,5 +575,5 @@ continue; } - Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass); + Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass, NULL); if (isnan(Mgrp)) { Ngrp ++; @@ -661,5 +661,5 @@ Mmos = getMmos (m, c); if (isnan(Mmos)) continue; - Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass); + Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass, NULL); if (isnan(Mgrp)) continue; Mrel = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); Index: trunk/Ohana/src/relphot/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageOps.c (revision 41461) +++ trunk/Ohana/src/relphot/src/ImageOps.c (revision 41462) @@ -581,10 +581,10 @@ /* determine Mcal values for all images */ -void setMcal (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { +void setMcal (Catalog *catalog, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr) { off_t i, j, m, c, n; int mark, bad, Nfew, Nbad, Nmos, Ngrp, Nrel, Ngrid, Nsys; - // FREEZE_IMAGES only applies to mosaic data (eg, gpc1) + if (IMAGE_ZPT_MODE == IMAGE_ZPT_MODE_NONE) return; fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE); @@ -609,4 +609,12 @@ FitDataSetAlloc (&brightStars, Nmax, 0, 0); + // until the analysis has converged a bit, do not use the IRLS analysis + if (Nloop < 3) { + brightStars.MaxIterations = 0; + kronStars.MaxIterations = 0; + psfStars.MaxIterations = 0; + } + + Nfew = Nbad = Nmos = Ngrp = Ngrid = Nrel = Nsys = 0; @@ -627,9 +635,26 @@ } - // only freeze the gpc1 chips - // XXX : when do we FREEZE_IMAGES in a mix of GPC1, GPC2, etc? - // if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue; - if (FREEZE_IMAGES) continue; + // if requested, freeze mosaic chips: + if (FREEZE_IMAGES && isMosaicChip (image[i].photcode)) continue; if (image[i].photcode == 0) continue; // skip the PHU images + + int badNight = (image[i].flags & ID_IMAGE_NIGHT_POOR); + int badMosaic = (image[i].flags & ID_IMAGE_MOSAIC_POOR); + + // in BAD_NIGHT mode, we fit ONLY images in bad nights + if ((IMAGE_ZPT_MODE == IMAGE_ZPT_MODE_BAD_NIGHT)) { + if (!badNight) continue; + } + + // 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) continue; + } + + // in BAD_MOSAIC mode, we fit bad mosaics ignoring night state + if ((IMAGE_ZPT_MODE == IMAGE_ZPT_MODE_BAD_MOSAIC)) { + if (badMosaic) continue; + } // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0: @@ -665,5 +690,5 @@ continue; } - float Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass); + float Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass, NULL); // ignore error for now? if (isnan(Mgrp)) { Ngrp ++; @@ -1027,5 +1052,5 @@ Mmos = getMmos (m, c); if (isnan(Mmos)) continue; - Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass); + Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass, NULL); if (isnan(Mgrp)) continue; Mgrid = getMgrid (m, c); Index: trunk/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- trunk/Ohana/src/relphot/src/MosaicOps.c (revision 41461) +++ trunk/Ohana/src/relphot/src/MosaicOps.c (revision 41462) @@ -10,5 +10,5 @@ // relationships between the mosaics and their associated images static off_t *MosaicN_Image; // number of images associated with the given mosaic -static off_t **MosaicToImage; // list of imagesa associated with the given mosaic +static off_t **MosaicToImage; // list of images associated with the given mosaic // mosaic index for given image : ImageToMosaic[ImageIndex] = MosaicIndex (ImageIndex : 0 < Nimage) @@ -219,4 +219,5 @@ mosaic[i].photcode = 0; mosaic[i].skipCal = FALSE; + mosaic[i].inTGroup = FALSE; // not (yet?) assigned to a TGroup memset (&mosaic[i].coords, 0, sizeof(Coords)); @@ -301,4 +302,7 @@ mosaic[j].dMcal = 0.0; mosaic[j].McalChiSq = 0.0; + + mosaic[j].nFitPhotom = 0; + mosaic[j].nValPhotom = 0; mosaic[j].dMsys = subset[i].flags; @@ -750,4 +754,41 @@ } +void markBadMosaic (int myMosaic) { + + off_t Nimage; + Image *image; + + // this mosaic has been identified as poor. + mosaic[myMosaic].flags |= ID_IMAGE_PHOTOM_POOR; + mosaic[myMosaic].McalPSF = 0.0; + mosaic[myMosaic].McalAPER = 0.0; + + image = getimages (&Nimage, NULL); + + // all images should be marked as coming from a bad mosaic + // these will be fitted independently + for (off_t i = 0; i < MosaicN_Image[myMosaic]; i++) { + off_t im = MosaicToImage[myMosaic][i]; + image[im].flags |= ID_IMAGE_MOSAIC_POOR; + } +} + +void markGoodMosaic (int myMosaic) { + + off_t Nimage; + Image *image; + + // this mosaic has been identified as good. + mosaic[myMosaic].flags &= ~ID_IMAGE_PHOTOM_POOR; + + image = getimages (&Nimage, NULL); + + // all images should be marked as NOT coming from a bad mosaic + for (off_t i = 0; i < MosaicN_Image[myMosaic]; i++) { + off_t im = MosaicToImage[myMosaic][i]; + image[im].flags &= ~ID_IMAGE_MOSAIC_POOR; + } +} + void initMosaicBins (Catalog *catalog, int Ncatalog, int doMosaicList) { @@ -923,4 +964,5 @@ int Nsys; int Nskip; + int Nloop; off_t Nmax; int PoorImages; @@ -943,7 +985,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, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr); + +void SetMmosInfoInit (SetMmosInfo *info, off_t Nmax, int allocLists, int PoorImages, int Nloop) { info->Nfew = 0; info->Nbad = 0; @@ -955,4 +997,5 @@ info->Nmax = Nmax; + info->Nloop = Nloop; info->PoorImages = PoorImages; @@ -962,4 +1005,11 @@ FitDataSetAlloc (&info->kronStars, Nmax, 0, 0); FitDataSetAlloc (&info->brightStars, Nmax, 0, 0); + + // until the analysis has converged a bit, do not use the IRLS analysis + if (Nloop < 3) { + info->brightStars.MaxIterations = 0; + info->kronStars.MaxIterations = 0; + info->psfStars.MaxIterations = 0; + } } } @@ -1003,16 +1053,18 @@ } -int setMmos (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { +int setMmos (Catalog *catalog, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr) { off_t i, N, Nmax; Image *image; - if (!MOSAIC_ZEROPT) return (FALSE); - if (FREEZE_MOSAICS) return (FALSE); + if (!MOSAIC_ZEROPT) return FALSE; + if (FREEZE_MOSAICS) return FALSE; + + if (MOSAIC_ZPT_MODE == MOSAIC_ZPT_MODE_NONE) return FALSE; // plots cannot be done in a threaded context (the plot commands collide) // 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, Nloop, PoorImages, flatcorr); return status; } @@ -1034,5 +1086,5 @@ SetMmosInfo info; - SetMmosInfoInit (&info, Nmax, TRUE, PoorImages); + SetMmosInfoInit (&info, Nmax, TRUE, PoorImages, Nloop); // int savePlotDelay = PLOTDELAY; @@ -1078,4 +1130,29 @@ assert (Nmos >= 0); assert (Nmos < Nmosaic); + + int badNight = (myMosaic[0].flags & ID_IMAGE_NIGHT_POOR); + int badMosaic = (myMosaic[0].flags & ID_IMAGE_MOSAIC_POOR); + + // 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 + if (MOSAIC_ZPT_MODE == MOSAIC_ZPT_MODE_BAD_NIGHT) { + if (!badNight) return TRUE; + } + + // 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; + } + + // in GOOD_MOSAIC mode, we fit good mosaics ignoring night state + if ((MOSAIC_ZPT_MODE == MOSAIC_ZPT_MODE_GOOD_MOSAIC)) { + if (badMosaic) return TRUE; + } // Image *imageReal = getimages (&NimageReal, NULL); returned pointer is not used @@ -1140,5 +1217,5 @@ if (fout) { float Mcal = getMcal (m, c, MAG_CLASS_PSF); - // float Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass); + // 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); @@ -1164,5 +1241,5 @@ continue; } - float Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass); + float Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass, NULL); if (isnan(Mgrp)) { info->Ngrp++; @@ -1257,4 +1334,5 @@ myMosaic[0].McalChiSq = psfStars->chisq; myMosaic[0].nFitPhotom = psfStars->Nmeas; + myMosaic[0].nValPhotom = N; // fit1d_irls (kronStars, N); // does it make sense to calculate a zero point based kron mags? @@ -1296,5 +1374,5 @@ } -int setMmos_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { +int setMmos_threaded (Catalog *catalog, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr) { int i; @@ -1317,5 +1395,5 @@ SetMmosInfo summary; - SetMmosInfoInit (&summary, Nmax, FALSE, PoorImages); + SetMmosInfoInit (&summary, Nmax, FALSE, PoorImages, Nloop); pthread_attr_t attr; @@ -1340,5 +1418,8 @@ threadinfo[i].image = image; threadinfo[i].flatcorr = flatcorr; - SetMmosInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages); + + // 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, Nloop); pthread_create (&threads[i], NULL, setMmos_worker, &threadinfo[i]); } @@ -1398,5 +1479,5 @@ SetMmosInfo results; - SetMmosInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages); // allocate list, dlist arrays here + SetMmosInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages, threadinfo->info.Nloop); while (1) { @@ -1709,5 +1790,5 @@ Mcal = getMcal (m, c, MAG_CLASS_PSF); if (isnan(Mcal)) continue; - Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass); + Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass, NULL); if (isnan(Mgrp)) continue; Mgrid = getMgrid (m, c); @@ -1767,5 +1848,5 @@ void clean_mosaics () { - off_t i, N, mark, Nmark, Nscatter, Noffset; + off_t i, N, mark, Ntotal, Nmark, Nscatter, Noffset, Nrejects, Nchisq; double *mlist, *slist, *dlist; double MaxOffset, MedOffset, MaxScatter; @@ -1800,10 +1881,17 @@ fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset); - Nmark = Nscatter = Noffset = 0; + // XXX this needs to be configurable + float MaxChiSq = 10.0; + + Ntotal = Nmark = Nscatter = Noffset = Nrejects = Nchisq = 0; for (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; - if (mosaic[i].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP)) continue; + Ntotal ++; + + if (mosaic[i].flags & (ID_IMAGE_PHOTOM_FEW)) { + continue; + } if (mosaic[i].skipCal) continue; @@ -1817,13 +1905,21 @@ Noffset ++; } + if (mosaic[i].nFitPhotom < 0.7*mosaic[i].nValPhotom) { + mark = TRUE; + Nrejects ++; + } + if (mosaic[i].McalChiSq > MaxChiSq) { + mark = TRUE; + Nchisq ++; + } if (mark) { Nmark ++; - mosaic[i].flags |= ID_IMAGE_PHOTOM_POOR; + markBadMosaic(i); // mark the images associated with a bad night } else { - mosaic[i].flags &= ~ID_IMAGE_PHOTOM_POOR; - } - } - - fprintf (stderr, OFF_T_FMT" mosaics marked poor ("OFF_T_FMT" scatter, "OFF_T_FMT" offset)\n", Nmark, Nscatter, Noffset); + markGoodMosaic(i); // mark the images associated with a good night + } + } + + 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); free (mlist); Index: trunk/Ohana/src/relphot/src/StarOps.c =================================================================== --- trunk/Ohana/src/relphot/src/StarOps.c (revision 41461) +++ trunk/Ohana/src/relphot/src/StarOps.c (revision 41462) @@ -493,5 +493,5 @@ float Mmos = getMmos (m, i); if (isnan(Mmos)) continue; - float Mgrp = getMgrp (m, i, catalog[i].measureT[m].airmass); + float Mgrp = getMgrp (m, i, catalog[i].measureT[m].airmass, NULL); if (isnan(Mgrp)) continue; float Mgrid = getMgrid (m, i); @@ -733,5 +733,5 @@ Mmos = getMmos (m, i); if (isnan(Mmos)) { Nmos ++; continue; } - Mgrp = getMgrp (m, i, catalog[i].measureT[m].airmass); + Mgrp = getMgrp (m, i, catalog[i].measureT[m].airmass, NULL); if (isnan(Mgrp)) { Ngrp ++; continue; } Mgrid = getMgrid (m, i); @@ -1082,2 +1082,37 @@ } +void dump_catalog (Catalog *catalog, off_t c, int Npass) { + + int Nsecfilt = GetPhotcodeNsecfilt (); + + FILE *fout = NULL; + char filename[64]; + snprintf (filename, 64, "tcat.%02d.dat", Npass); + fout = fopen (filename, "w"); + + for (off_t n = 0; n < catalog[c].Naverage; n++) { + + off_t m = catalog[c].averageT[n].measureOffset; + for (off_t k = 0; k < catalog[c].averageT[n].Nmeasure; k++, m++) { + + float Mcal = getMcal (m, c, MAG_CLASS_PSF); // image zero point + float Mmos = getMmos (m, c); // mosaic zero point + float Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass, NULL); + + float MrelPSF = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude + // float Mgrid = getMgrid (m, c); // camera offset (deprecated?) + // float Mflat = getMflat (m, c, flatcorr, catalog); // flat-field correction + + off_t n = catalog[c].measureT[m].averef; + + // nominal magnitude for this measurement (instrumental + Ko(airmass - 1) + Co + 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 - Mgrp; + + fprintf (fout, "%f %f : %f %f : %f %f %f : %f %f\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, Mcal, Mmos, Mgrp, catalog[c].measureT[m].airmass, delta); + } + } + fclose (fout); +} Index: trunk/Ohana/src/relphot/src/TGroupOps.c =================================================================== --- trunk/Ohana/src/relphot/src/TGroupOps.c (revision 41461) +++ trunk/Ohana/src/relphot/src/TGroupOps.c (revision 41462) @@ -1,3 +1,5 @@ # include "relphot.h" + +void assignMosaicsToTGroups (void); // tgroupTimes carries the times of the tgroups (initially, just the photometric nights) @@ -78,4 +80,9 @@ tgroup[j].dMsys = 0.0; tgroup[j].McalChiSq = 0.0; // NAN or 0.0? + + tgroup[j].stdev = NAN; + tgroup[j].nFitPhotom = 0; + tgroup[j].nValPhotom = 0; + tgroup[j].flags = ID_IMAGE_PHOTOM_UBERCAL; // we start by treating the suspected photometric nights as photometric tgroup[j].photcode = photcodes[j][0].code; @@ -85,4 +92,8 @@ ALLOCATE (tgroup[j].image, off_t, tgroup[j].NIMAGE); tgroup[j].image[0] = -1; + + tgroup[j].NMOSAIC = 100; + tgroup[j].Nmosaic = 0; + ALLOCATE (tgroup[j].mosaic, Mosaic *, tgroup[j].NMOSAIC); tgroup[j].Nmeasure = 0; @@ -129,10 +140,21 @@ tgroup[j].dMsys = 0.0; tgroup[j].McalChiSq = 0.0; // NAN or 0.0? + + tgroup[j].stdev = NAN; + tgroup[j].nFitPhotom = 0; + tgroup[j].nValPhotom = 0; + tgroup[j].flags = 0; // tgroups NOT in the original file are not suspected as photometric tgroup[j].photcode = photcodes[j][0].code; + tgroup[j].NIMAGE = 1000; tgroup[j].Nimage = 0; // the new TGroup is not suspected to be photometric ALLOCATE (tgroup[j].image, off_t, tgroup[j].NIMAGE); tgroup[j].image[0] = -1; + + tgroup[j].NMOSAIC = 100; + tgroup[j].Nmosaic = 0; + ALLOCATE (tgroup[j].mosaic, Mosaic *, tgroup[j].NMOSAIC); + tgroup[j].Nmeasure = 0; tgroup[j].NMEASURE = 1000; @@ -199,4 +221,7 @@ } + // assign the mosaics to each tgroup -- we do this by image + assignMosaicsToTGroups (); + initTGroupsMcal (); @@ -215,4 +240,5 @@ for (int j = 0; j < tgroupTimes[i][0].nCode; j++) { free (tgroup[j].image); + free (tgroup[j].mosaic); free (tgroup[j].measure); free (tgroup[j].catalog); @@ -249,4 +275,8 @@ tgroup[j].dMsys = 0.0; tgroup[j].McalChiSq = 0.0;// NAN or 0.0? + + tgroup[j].stdev = NAN; + tgroup[j].nFitPhotom = 0; + tgroup[j].nValPhotom = 0; if (tgroup[j].Nimage == 0) continue; // no images, ignore the tgroup @@ -278,4 +308,38 @@ } +void assignMosaicsToTGroups (void) { + + if (!TGROUP_ZEROPT) return; + + // init the tgroup array values + 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].Nimage == 0) continue; // no images, ignore the tgroup + + // find the mosaics by image (if MOSAIC_ZEROPT is active) + for (int im = 0; MOSAIC_ZEROPT && (im < tgroup[j].Nimage); im++) { + int imageIdx = tgroup[j].image[im]; + Mosaic *mosaic = getMosaicForImage (imageIdx); + + if (mosaic->inTGroup) continue; + + // mosaic to tgroup and extend as needed + tgroup[j].mosaic[tgroup[j].Nmosaic] = mosaic; + mosaic->inTGroup = TRUE; + + // advance array and extend if needed + tgroup[j].Nmosaic ++; + if (tgroup[j].Nmosaic >= tgroup[j].NMOSAIC) { + tgroup[j].NMOSAIC += 100; + REALLOCATE (tgroup[j].mosaic, Mosaic *, tgroup[j].NMOSAIC); + } + } + } + } + return; +} + TGroup *getTGroupForImage (off_t im) { @@ -345,4 +409,6 @@ fprintf (stderr, "*** return Mcal from tgroup.Mcal to image.Mcal ***\n"); + + // NOTE: nights which are poor should have had Mcal set to 0, making this a NOP for those nights // copy the tgroup results to the images. set the tgroup Mcal to 0.0 since we have moved its @@ -360,4 +426,6 @@ image[im].dMagSys = tgroup[i].dMsys; image[im].nFitPhotom = tgroup[i].nFitPhotom; + image[im].flags |= (tgroup[i].flags & ID_IMAGE_NIGHT_POOR); // XXX copy this down? + image[im].flags |= (tgroup[i].flags & ID_IMAGE_MOSAIC_POOR); // XXX copy this down? image[im].flags |= (tgroup[i].flags & ID_IMAGE_PHOTOM_FEW); image[im].flags |= (tgroup[i].flags & ID_IMAGE_PHOTOM_POOR); @@ -371,4 +439,52 @@ } } +} + +void markBadTGroup (TGroup *tgroup) { + + off_t Nimage; + Image *image; + + // this tgroup has been identified as poor. + tgroup->flags |= ID_IMAGE_NIGHT_POOR; + + // all associated mosaics should be marked as coming from a bad night + // these will be fitted independently + for (off_t j = 0; j < tgroup->Nmosaic; j++) { + tgroup->mosaic[j]->flags |= ID_IMAGE_NIGHT_POOR; + } + + // all associated images should be marked as coming from a bad night + // these will be fitted independently + image = getimages (&Nimage, NULL); + for (off_t j = 0; j < tgroup->Nimage; j++) { + off_t im = tgroup->image[j]; + image[im].flags |= ID_IMAGE_NIGHT_POOR; + } + + tgroup->McalPSF = 0.0; + tgroup->McalAPER = 0.0; + tgroup->dKlam = 0.0; +} + +void markGoodTGroup (TGroup *tgroup) { + + off_t Nimage; + Image *image; + + // this tgroup has been identified as good. + tgroup->flags &= ~ID_IMAGE_NIGHT_POOR; + + // all associated mosaics should be marked as NOT coming from a bad night + for (off_t j = 0; j < tgroup->Nmosaic; j++) { + tgroup->mosaic[j]->flags &= ~ID_IMAGE_NIGHT_POOR; + } + + // all associated images should be marked as NOT coming from a bad night + image = getimages (&Nimage, NULL); + for (off_t j = 0; j < tgroup->Nimage; j++) { + off_t im = tgroup->image[j]; + image[im].flags &= ~ID_IMAGE_NIGHT_POOR; + } } @@ -488,17 +604,18 @@ } -float getMgrp (off_t meas, int cat, float airmass) { - - float value; - - if (!TGROUP_ZEROPT) return (0); +// XXX return the error as well as the value +float getMgrp (off_t meas, int cat, float airmass, float *dZpt) { + + if (dZpt) { *dZpt = 0; } + + if (!TGROUP_ZEROPT) return 0.0; // unassigned measurements belong to simple chips TGroup *myGroup = MeasureToTGroup[cat][meas]; - if (!myGroup) return (0.0); - - // if (myGroup->flags & IMAGE_BAD) return (NAN); - value = myGroup->McalPSF + myGroup->dKlam*(airmass - 1.0); - return (value); + if (!myGroup) return 0.0; + + float value = myGroup->McalPSF + myGroup->dKlam*(airmass - 1.0); + if (dZpt) { *dZpt = myGroup->dMcal; } + return value; } @@ -525,4 +642,5 @@ int Nsys; int Nskip; + int Nloop; off_t Nmax; int PoorImages; @@ -545,7 +663,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, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr); + +void SetMgrpInfoInit (SetMgrpInfo *info, off_t Nmax, int allocLists, int PoorImages, int Nloop) { info->Nfew = 0; info->Nbad = 0; @@ -558,4 +676,5 @@ info->Nmax = Nmax; + info->Nloop = Nloop; info->PoorImages = PoorImages; @@ -567,4 +686,10 @@ FitDataSetAlloc (&info->brightStars, Nmax, fitOrder, 0); + // until the analysis has converged a bit, do not use the IRLS analysis + if (Nloop < 100) { + info->brightStars.MaxIterations = 0; + info->kronStars.MaxIterations = 0; + info->psfStars.MaxIterations = 0; + } if (1) { FitDataSetAddPriors (&info->psfStars); @@ -614,5 +739,5 @@ } -int setMgrp (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { +int setMgrp (Catalog *catalog, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr) { off_t N; @@ -620,19 +745,14 @@ if (!TGROUP_ZEROPT) return (FALSE); + if (TGROUP_ZPT_MODE == TGROUP_ZPT_MODE_NONE) return (FALSE); // plots cannot be done in a threaded context (the plot commands collide) // 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, Nloop, PoorImages, flatcorr); return status; } image = getimages (&N, NULL); - - 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; @@ -645,8 +765,5 @@ SetMgrpInfo info; - SetMgrpInfoInit (&info, Nmax, TRUE, PoorImages); - - // int savePlotDelay = PLOTDELAY; - // if (PLOTSTUFF) PLOTDELAY = 0.0; + SetMgrpInfoInit (&info, Nmax, TRUE, PoorImages, Nloop); for (off_t i = 0; i < NtgroupTimes; i++) { @@ -654,5 +771,5 @@ 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, 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, 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); } } @@ -668,9 +785,4 @@ fprintf (stderr, "%d tgroups 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); } @@ -684,9 +796,10 @@ off_t j, NimageReal; - /* on PoorImages run, skip good images */ - if (info->PoorImages) { - int bad = myTGroup->flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP); - if (!bad) return TRUE; - } + // 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; + } FitDataSet *psfStars = &info->psfStars; @@ -704,5 +817,6 @@ int testImage = FALSE; TGTimes *tgroup = (TGTimes *) myTGroup->parent; - testImage |= (abs(tgroup->start - 1245283200) < 10); + // testImage |= (abs(tgroup->start - 1245283200) < 10); + // testImage = TRUE; FILE *fout = NULL; @@ -825,19 +939,17 @@ /* too few good measurements or too many bad measurements (skip in PoorImages run) */ - if (!info->PoorImages) { - int mark = (N < NIGHT_TOOFEW) || (N < NIGHT_GOOD_FRACTION*myTGroup->Nmeasure); - if (mark) { - if (VERBOSE2) { - 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); - } - myTGroup->flags |= ID_IMAGE_PHOTOM_FEW; - info->Nfew ++; - if (testImage) { - fprintf (stderr, "NOTE: *** marked test image poor : %d %d %d***\n", (int) N, (int) NIGHT_TOOFEW, (int) (NIGHT_GOOD_FRACTION*myTGroup->Nmeasure)); - } - } else { - myTGroup->flags &= ~ID_IMAGE_PHOTOM_FEW; - } + int mark = (N < NIGHT_TOOFEW) || (N < NIGHT_GOOD_FRACTION*myTGroup->Nmeasure); + if (mark) { + if (VERBOSE2) { + 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); + } + myTGroup->flags |= ID_IMAGE_PHOTOM_FEW; + info->Nfew ++; + if (testImage) { + fprintf (stderr, "NOTE: *** marked test image poor : %d %d %d***\n", (int) N, (int) NIGHT_TOOFEW, (int) (NIGHT_GOOD_FRACTION*myTGroup->Nmeasure)); + } + } else { + myTGroup->flags &= ~ID_IMAGE_PHOTOM_FEW; } @@ -852,6 +964,8 @@ myTGroup->McalPSF = psfStars->bSaveArray[0][0]; myTGroup->dMcal = psfStars->bSigma[0]; + myTGroup->stdev = psfStars->sigma; myTGroup->McalChiSq = psfStars->chisq; myTGroup->nFitPhotom = psfStars->Nmeas; + myTGroup->nValPhotom = N; myTGroup->dKlam = psfStars->bSaveArray[1][0]; @@ -864,4 +978,6 @@ fprintf (stderr, "test night %.0f aper: %f %f %d ... ", ohana_sec_to_mjd(parent->start), myTGroup->McalPSF, myTGroup->dMcal, myTGroup->nFitPhotom); } + + // fprintf (stderr, "Mgrp: %6.3f +/- %6.3f %5d of %5d | %.2f %.0f\n", myTGroup->McalPSF, myTGroup->dMcal, myTGroup->nFitPhotom, myTGroup->nValPhotom, myTGroup->McalChiSq, ohana_sec_to_mjd(tgroup->start)); if (PLOTSTUFF) { @@ -885,12 +1001,8 @@ } - // FitDataSetFree (&brightStars); - // FitDataSetFree (&kronStars); - // FitDataSetFree (&psfStars); - return TRUE; } -int setMgrp_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { +int setMgrp_threaded (Catalog *catalog, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr) { int i; @@ -901,10 +1013,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; for (off_t i = 0; i < NtgroupTimes; i++) { @@ -916,5 +1022,5 @@ SetMgrpInfo summary; - SetMgrpInfoInit (&summary, Nmax, FALSE, PoorImages); + SetMgrpInfoInit (&summary, Nmax, FALSE, PoorImages, Nloop); pthread_attr_t attr; @@ -939,5 +1045,8 @@ threadinfo[i].image = image; threadinfo[i].flatcorr = flatcorr; - SetMgrpInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages); + + // 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 + SetMgrpInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages, Nloop); pthread_create (&threads[i], NULL, setMgrp_worker, &threadinfo[i]); } @@ -986,9 +1095,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; } @@ -999,5 +1103,5 @@ SetMgrpInfo results; - SetMgrpInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages); // allocate list, dlist arrays here + SetMgrpInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages, threadinfo->info.Nloop); // allocate list, dlist arrays here while (1) { @@ -1024,10 +1128,9 @@ } -# if (0) +// XXX should I split these stats out by photcode? +// XXX double check the IMAGE_BAD filter: StatType statsTGroupM (Catalog *catalog) { OHANA_UNUSED_PARAM(catalog); - off_t i; - int n; double *list, *dlist; StatType stats; @@ -1037,12 +1140,15 @@ // XXX do this by time & photcode? - ALLOCATE (list, double, Ntgroup); - ALLOCATE (dlist, double, Ntgroup); - - n = 0; + ALLOCATE (list, double, NtgroupTimes*Nphotcodes); + 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 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 & NIGHT_BAD) continue; list[n] = tgroup[j].McalPSF; dlist[n] = 1; @@ -1062,5 +1168,4 @@ OHANA_UNUSED_PARAM(catalog); - off_t i, n; double *list, *dlist; StatType stats; @@ -1069,13 +1174,16 @@ if (!TGROUP_ZEROPT) return (stats); - ALLOCATE (list, double, Ntgroup); - ALLOCATE (dlist, double, Ntgroup); - - n = 0; - for (i = 0; i < Ntgroup; i++) { - if (tgroup[i].flags & IMAGE_BAD) continue; - list[n] = tgroup[i].dMcal; - dlist[n] = 1; - n++; + ALLOCATE (list, double, NtgroupTimes*Nphotcodes); + ALLOCATE (dlist, double, NtgroupTimes*Nphotcodes); + + 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; + list[n] = tgroup[j].dMcal; + dlist[n] = 1; + n++; + } } @@ -1088,9 +1196,8 @@ } -StatType statsTGroupN (Catalog *catalog, FlatCorrectionTable *flatcorr) { - - off_t i, j, m, c, n, N; +StatType statsTGroupX (Catalog *catalog) { + OHANA_UNUSED_PARAM(catalog); + double *list, *dlist; - float Mcal, Mgrid, Mrel, Mmos; StatType stats; @@ -1098,32 +1205,17 @@ if (!TGROUP_ZEROPT) return (stats); - ALLOCATE (list, double, Ntgroup); - ALLOCATE (dlist, double, Ntgroup); - - n = 0; - for (i = 0; i < Ntgroup; i++) { - if (tgroup[i].flags & IMAGE_BAD) continue; - - N = 0; - for (j = 0; j < tgroup[i].Nmeasure; j++) { - - m = tgroup[i].measure[j]; - c = tgroup[i].catalog[j]; - - Mcal = getMcal (m, c, MAG_CLASS_PSF); - if (isnan(Mcal)) continue; - Mmos = getMmos (m, c); - if (isnan(Mmos)) continue; - Mgrid = getMgrid (m, c); - if (isnan(Mgrid)) continue; - Mrel = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); - if (isnan(Mrel)) continue; - N++; - } - list[n] = N; - dlist[n] = 1; - n++; - } - // fprintf (stderr, "Ntgroup: "OFF_T_FMT", n: "OFF_T_FMT"\n", Ntgroup, n); + ALLOCATE (list, double, NtgroupTimes*Nphotcodes); + ALLOCATE (dlist, double, NtgroupTimes*Nphotcodes); + + 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; + list[n] = tgroup[j].McalChiSq; + dlist[n] = 1; + n++; + } + } liststats_setmode (&stats, "MEAN"); @@ -1135,40 +1227,8 @@ } -StatType statsTGroupX (Catalog *catalog) { - OHANA_UNUSED_PARAM(catalog); - - off_t i, n; - double *list, *dlist; - StatType stats; - - bzero (&stats, sizeof (StatType)); - if (!TGROUP_ZEROPT) return (stats); - - ALLOCATE (list, double, Ntgroup); - ALLOCATE (dlist, double, Ntgroup); - - n = 0; - for (i = 0; i < Ntgroup; i++) { - if (tgroup[i].flags & IMAGE_BAD) continue; - if (KEEP_UBERCAL && (tgroup[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; - - list[n] = tgroup[i].McalChiSq; - dlist[n] = 1; - n++; - } - - liststats_setmode (&stats, "MEAN"); - - liststats (list, dlist, NULL, n, &stats); - free (list); - free (dlist); - return (stats); -} -# endif - /* mark tgroup if: abs(Mcal - ) too large, dMcal too large */ void clean_tgroups () { - off_t i, N, mark, Nmark, Nscatter, Noffset; + off_t mark, Ntotal, Nmark, Nscatter, Noffset, Nrejects, Nchisq; double *mlist, *slist, *dlist; double MaxOffset, MedOffset, MaxScatter; @@ -1182,5 +1242,6 @@ ALLOCATE (dlist, double, NtgroupTimes*Nphotcodes); - for (i = N = 0; i < NtgroupTimes; i++) { + 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++) { @@ -1197,40 +1258,56 @@ liststats (mlist, dlist, NULL, N, &stats); - MaxOffset = MAX (NIGHT_OFFSET, 2*stats.sigma); + 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, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset); + fprintf (stderr, "TGROUPS: median Mrel: %f, median dMrel: %f, Max Offset: %f, Max Scatter: %f\n", MedOffset, stats.median, MaxOffset, MaxScatter); - Nmark = Nscatter = Noffset = 0; - for (i = 0; i < NtgroupTimes; i++) { + // XXX this needs to be configurable + float MaxChiSq = 10.0; + + Ntotal = Nmark = Nscatter = Noffset = Nrejects = Nchisq = 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 we are keeping ubercal sacrosanct, then we should not be allowed to break them... - if (KEEP_UBERCAL && (tgroup[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; - if (tgroup[i].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP)) continue; + 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 + if (tgroup[j].flags & ID_IMAGE_PHOTOM_FEW) { + markBadTGroup(&tgroup[j]); // mark the image associated with a bad night + continue; + } mark = FALSE; - if (tgroup[i].dMcal > MaxScatter) { + if (tgroup[j].dMcal > MaxScatter) { mark = TRUE; Nscatter ++; } - if (fabs(tgroup[i].McalPSF - MedOffset) > MaxOffset) { + if (fabs(tgroup[j].McalPSF - MedOffset) > MaxOffset) { mark = TRUE; Noffset ++; } + if (tgroup[j].nFitPhotom < 0.7*tgroup[j].nValPhotom) { + mark = TRUE; + Nrejects ++; + } + if (tgroup[j].McalChiSq > MaxChiSq) { + mark = TRUE; + Nchisq ++; + } if (mark) { Nmark ++; - tgroup[i].flags |= ID_IMAGE_PHOTOM_POOR; + markBadTGroup(&tgroup[j]); // mark the tgroup & associated images & mosaics with a bad night } else { - tgroup[i].flags &= ~ID_IMAGE_PHOTOM_POOR; + markGoodTGroup(&tgroup[j]); // mark the tgroup & associated images & mosaics with a good night } - // XXX : unset the ubercal flag if we mark the night as bad? - } - } - - fprintf (stderr, OFF_T_FMT" tgroups marked poor ("OFF_T_FMT" scatter, "OFF_T_FMT" offset)\n", Nmark, Nscatter, Noffset); + } + } + + 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); free (mlist); @@ -1347,2 +1424,47 @@ } +void dump_tgroups (Catalog *catalog, int Npass) { + + int Nsecfilt = GetPhotcodeNsecfilt (); + + int Ngrp = 0; + 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++) { + + FILE *fout = NULL; + char filename[64]; + snprintf (filename, 64, "tgrp.%05d.%02d.dat", Ngrp, Npass); + fout = fopen (filename, "w"); + Ngrp ++; + + TGroup *myTGroup = &tgroup[j]; + + for (int k = 0; k < myTGroup->Nmeasure; k++) { + + off_t m = myTGroup->measure[k]; + off_t c = myTGroup->catalog[k]; + + float Mcal = getMcal (m, c, MAG_CLASS_PSF); // image zero point + float Mmos = getMmos (m, c); // mosaic zero point + float Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass, NULL); + + float MrelPSF = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude + // float Mgrid = getMgrid (m, c); // camera offset (deprecated?) + // float Mflat = getMflat (m, c, flatcorr, catalog); // flat-field correction + + off_t n = catalog[c].measureT[m].averef; + + // nominal magnitude for this measurement (instrumental + Ko(airmass - 1) + Co + 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 - Mgrp; + + fprintf (fout, "%f %f : %f %f : %f %f %f : %f %f\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, Mcal, Mmos, Mgrp, catalog[c].measureT[m].airmass, delta); + } + fclose (fout); + } + } +} Index: trunk/Ohana/src/relphot/src/ZeroPointModes.c =================================================================== --- trunk/Ohana/src/relphot/src/ZeroPointModes.c (revision 41462) +++ trunk/Ohana/src/relphot/src/ZeroPointModes.c (revision 41462) @@ -0,0 +1,75 @@ +# include "relphot.h" + +void SetZeroPointModes (int nloop) { + + if (TGROUP_ZEROPT && MOSAIC_ZEROPT) { + if ((nloop > 4) && (nloop <= 8)) { + // after iterating a few times on the TGroup zero points: + // * identify the photometric nights + // * fit the mosaics from the non-photometric nights + TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_GOOD_NIGHT; // stop fitting the bad nights + MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_BAD_NIGHT; + IMAGE_ZPT_MODE = IMAGE_ZPT_MODE_NONE; + clean_tgroups(); // do this on every pass to update the status of nights + } + if ((nloop > 8) && (nloop <= 12)) { + // 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 + } + } + + if (TGROUP_ZEROPT && !MOSAIC_ZEROPT) { + if ((nloop > 4) && (nloop <= 8)) { + // after iterating a few times on the TGroup zero points: + // * identify the photometric nights + // * fit the images from the non-photometric nights + TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_GOOD_NIGHT; // stop fitting the bad nights + MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE; + IMAGE_ZPT_MODE = IMAGE_ZPT_MODE_BAD_NIGHT; + clean_tgroups(); // do this on every pass or just sometimes? + } + } + + if (!TGROUP_ZEROPT && MOSAIC_ZEROPT) { + if ((nloop > 4) && (nloop <= 8)) { + // after iterating a few times on the Mosaic zero points: + // * identify the good / bad mosaics + // * fit the images from the bad mosaics + TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_NONE; // stop fitting the bad nights + MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_GOOD_MOSAIC; + IMAGE_ZPT_MODE = IMAGE_ZPT_MODE_BAD_MOSAIC; + clean_mosaics(); // do this on every pass or just sometimes? + } + } + + if (!TGROUP_ZEROPT && !MOSAIC_ZEROPT) { + if ((nloop > 4) && (nloop <= 8)) { + // 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 (i < NLOOP - 1) rationalize_mosaics (catalog, Ncatalog); + // if (i % 6 == 1) rationalize_images (); + + // NOTE : in the past, I was not iterating enough before cleaning. make sure we do + // at least 8 loops first -- that should get the systematic errors down to the ~1% + // level, even in cases where we have an even split between photometric data and + // data with 1 mag of extinction. + + // if ((i > 8) && (i % 8 == 2)) clean_measures (catalog, Ncatalog, FALSE, flatcorr); + // if ((i > 8) && (i % 8 == 3)) clean_stars (catalog, Ncatalog); + // if ((i > 8) && (i % 8 == 5)) clean_mosaics (); + // if ((i > 8) && (i % 8 == 5)) clean_tgroups (); + // if ((i > 8) && (i % 8 == 5)) clean_images (); + Index: trunk/Ohana/src/relphot/src/args.c =================================================================== --- trunk/Ohana/src/relphot/src/args.c (revision 41461) +++ trunk/Ohana/src/relphot/src/args.c (revision 41462) @@ -303,4 +303,33 @@ } + // global restriction on fitting image zero points for mosaic cameras + FREEZE_IMAGES = FALSE; + if ((N = get_argument (argc, argv, "-imfreeze"))) { + remove_argument (N, &argc, argv); + FREEZE_IMAGES = TRUE; + } + // global restriction on fitting mosaic zero points + FREEZE_MOSAICS = FALSE; + if ((N = get_argument (argc, argv, "-mosfreeze"))) { + remove_argument (N, &argc, argv); + FREEZE_MOSAICS = TRUE; + } + CALIBRATE_STACKS_AND_WARPS = FALSE; + if ((N = get_argument (argc, argv, "-only-stacks-and-warps"))) { + remove_argument (N, &argc, argv); + CALIBRATE_STACKS_AND_WARPS = TRUE; + } + USE_MCAL_PSF_FOR_STACK_APER = FALSE; + if ((N = get_argument (argc, argv, "-use-mcal-psf-for-stack-aper"))) { + remove_argument (N, &argc, argv); + USE_MCAL_PSF_FOR_STACK_APER = TRUE; + } + + + // default behavior is to fit all chips independently + TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_NONE; + MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE; + IMAGE_ZPT_MODE = IMAGE_ZPT_MODE_ALL; + /* XXX this argument used to do two things: specify the camera name and tell the analysis to calculate a common zero point for a single mosaic. I've moved the MOSAICNAME concept into the @@ -315,25 +344,6 @@ exit (2); } - } - - FREEZE_IMAGES = FALSE; - if ((N = get_argument (argc, argv, "-imfreeze"))) { - remove_argument (N, &argc, argv); - FREEZE_IMAGES = TRUE; - } - FREEZE_MOSAICS = FALSE; - if ((N = get_argument (argc, argv, "-mosfreeze"))) { - remove_argument (N, &argc, argv); - FREEZE_MOSAICS = TRUE; - } - CALIBRATE_STACKS_AND_WARPS = FALSE; - if ((N = get_argument (argc, argv, "-only-stacks-and-warps"))) { - remove_argument (N, &argc, argv); - CALIBRATE_STACKS_AND_WARPS = TRUE; - } - USE_MCAL_PSF_FOR_STACK_APER = FALSE; - if ((N = get_argument (argc, argv, "-use-mcal-psf-for-stack-aper"))) { - remove_argument (N, &argc, argv); - USE_MCAL_PSF_FOR_STACK_APER = TRUE; + IMAGE_ZPT_MODE = IMAGE_ZPT_MODE_NONE; + MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_ALL; } @@ -341,21 +351,17 @@ TGROUP_ZEROPT = FALSE; char *TGROUP_FILENAME = NULL; - if ((N = get_argument (argc, argv, "-tgroup"))) { + if ((N = get_argument (argc, argv, "-tgroup")) || (N = get_argument (argc, argv, "-tgroups"))) { TGROUP_ZEROPT = TRUE; remove_argument (N, &argc, argv); TGROUP_FILENAME = strcreate(argv[N]); remove_argument (N, &argc, argv); - FREEZE_IMAGES = TRUE; - } - if ((N = get_argument (argc, argv, "-tgroups"))) { - if (TGROUP_ZEROPT) { - fprintf (stderr, "use only one of -tgroups and -tgroup (same meaning)\n"); - exit (2); - } - TGROUP_ZEROPT = TRUE; - remove_argument (N, &argc, argv); - TGROUP_FILENAME = strcreate(argv[N]); - remove_argument (N, &argc, argv); - FREEZE_IMAGES = TRUE; + // if we fit for TGROUPS, start with the images and mosaics frozen + TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_ALL; + IMAGE_ZPT_MODE = IMAGE_ZPT_MODE_NONE; + MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE; + } + if ((N = get_argument (argc, argv, "-tgroup")) || (N = get_argument (argc, argv, "-tgroups"))) { + fprintf (stderr, "use only one of -tgroups and -tgroup (same meaning)\n"); + exit (2); } TGROUP_FIT_AIRMASS = FALSE; Index: trunk/Ohana/src/relphot/src/extra.c =================================================================== --- trunk/Ohana/src/relphot/src/extra.c (revision 41461) +++ trunk/Ohana/src/relphot/src/extra.c (revision 41462) @@ -130,2 +130,20 @@ } +double weight_cauchy (double x) { + double r = x / 2.385; + return (1.0 / (1.0 + SQ(r))); +} + +double VectorFractionInterpolate (double *values, float fraction, int Npts) { + + float F = fraction * Npts; + int N = fraction * Npts; + + if (N < 0 ) return NAN; + if (N >= Npts - 2) return NAN; + + // interpolate between N,N+1 + + double S = (F - N) * (values[N+1] - values[N]) + values[N]; + return S; +} Index: trunk/Ohana/src/relphot/src/fit1d_irls.c =================================================================== --- trunk/Ohana/src/relphot/src/fit1d_irls.c (revision 41461) +++ trunk/Ohana/src/relphot/src/fit1d_irls.c (revision 41462) @@ -2,5 +2,4 @@ // These should probably be tunable: -# define MAX_ITERATIONS 10 # define FIT_TOLERANCE 1e-4 # define FLT_TOLERANCE 1e-6 @@ -35,5 +34,5 @@ int converged = FALSE; - for (int iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations++) { + for (int iterations = 0; !converged && (iterations < dataset->MaxIterations); iterations++) { for (int i = 0; i < Npoints; i++) { @@ -293,4 +292,5 @@ dataset->mterm = 2*order + 1; dataset->Nbootstrap = Nbootstrap; + dataset->MaxIterations = 10; // arrays for actual data values Index: trunk/Ohana/src/relphot/src/global_stats.c =================================================================== --- trunk/Ohana/src/relphot/src/global_stats.c (revision 41461) +++ trunk/Ohana/src/relphot/src/global_stats.c (revision 41462) @@ -3,13 +3,15 @@ void global_stats (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr, int nloop) { - StatType stN, stX, stS, imN, imX, imM, imD, msM, msX, msN, msD; + StatType stN, stX, stS, imN, imX, imM, imD, msM, msX, msN, msD, tgM, tgX, tgD; liststats_init (&stN); liststats_init (&stX); liststats_init (&stS); + liststats_init (&imN); liststats_init (&imX); liststats_init (&imM); liststats_init (&imD); + liststats_init (&msN); liststats_init (&msX); @@ -17,5 +19,7 @@ liststats_init (&msD); - // INITTIME; + liststats_init (&tgX); + liststats_init (&tgM); + liststats_init (&tgD); fprintf (stderr, "\n"); @@ -29,11 +33,7 @@ int seccode = photcodes[Ns][0].code; - // MARKTIME("start star stats: %f sec\n", dtime); stN = statsStarN (catalog, Ncatalog, Nsec, seccode, flatcorr); - // MARKTIME("done star N: %f sec\n", dtime); stX = statsStarX (catalog, Ncatalog, Nsec); - // MARKTIME("done star X: %f sec\n", dtime); stS = statsStarS (catalog, Ncatalog, Nsec); - // MARKTIME("done star S: %f sec\n", dtime); fprintf (stderr, " --- stats for %s ---\n", photcodes[Ns][0].name); @@ -53,4 +53,8 @@ msX = statsMosaicX (catalog); + tgM = statsMosaicM (catalog); + tgD = statsMosaicdM (catalog); + tgX = statsMosaicX (catalog); + fprintf (stderr, "meas / image: %7.0f %7.0f %7.0f %7.0f %7.0f %6d\n", imN.median, imN.mean, imN.sigma, imN.min, imN.max, imN.Nmeas); fprintf (stderr, "Mcal image: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", imM.median, imM.mean, imM.sigma, imM.min, imM.max, imM.Nmeas); @@ -62,4 +66,8 @@ fprintf (stderr, "dMcal mosaic: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", msD.median, msD.mean, msD.sigma, msD.min, msD.max, msD.Nmeas); fprintf (stderr, "chisq mosaic: %7.1f %7.1f %7.1f %7.1f %7.1f %6d\n", msX.median, msX.mean, msX.sigma, msX.min, msX.max, msX.Nmeas); + + fprintf (stderr, "Mcal tgroup: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", tgM.median, tgM.mean, tgM.sigma, tgM.min, tgM.max, tgM.Nmeas); + fprintf (stderr, "dMcal tgroup: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", tgD.median, tgD.mean, tgD.sigma, tgD.min, tgD.max, tgD.Nmeas); + fprintf (stderr, "chisq tgroup: %7.1f %7.1f %7.1f %7.1f %7.1f %6d\n", tgX.median, tgX.mean, tgX.sigma, tgX.min, tgX.max, tgX.Nmeas); } Index: trunk/Ohana/src/relphot/src/liststats.c =================================================================== --- trunk/Ohana/src/relphot/src/liststats.c (revision 41461) +++ trunk/Ohana/src/relphot/src/liststats.c (revision 41462) @@ -417,23 +417,4 @@ } -double weight_cauchy (double x) { - double r = x / 2.385; - return (1.0 / (1.0 + SQ(r))); -} - -double VectorFractionInterpolate (double *values, float fraction, int Npts) { - - float F = fraction * Npts; - int N = fraction * Npts; - - if (N < 0 ) return NAN; - if (N >= Npts - 2) return NAN; - - // interpolate between N,N+1 - - double S = (F - N) * (values[N+1] - values[N]) + values[N]; - return S; -} - // we could define the weight to be the only scale factor: // \mu = \sum (value_i * weight_i) / \sum (weight_i) Index: trunk/Ohana/src/relphot/src/relphot_images.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_images.c (revision 41461) +++ trunk/Ohana/src/relphot/src/relphot_images.c (revision 41462) @@ -116,23 +116,16 @@ /* determine fit values */ for (i = 0; i < NLOOP; i++) { - if (PLOTSTUFF) { - plot_scatter (catalog, Ncatalog, flatcorr); - } - if (VERBOSE2) dumpMags (stderr, catalog, Ncatalog); - // checkImages("start loop"); + dump_tgroups (catalog, i); + dump_catalog (catalog, 0, i); // for a test, just dump a specific catalog setMrel (catalog, Ncatalog, flatcorr); // threaded - if (VERBOSE2) dumpMags (stderr, catalog, Ncatalog); - // checkImages("set Mrel"); - - if (PLOTSTUFF) { - plot_scatter (catalog, Ncatalog, flatcorr); - } - setMcal (catalog, FALSE, flatcorr); - // checkImages("set Mcal"); - if (VERBOSE2) dumpMags (stderr, catalog, Ncatalog); - - setMmos (catalog, FALSE, flatcorr); - setMgrp (catalog, FALSE, flatcorr); + dump_tgroups (catalog, i + 20); + dump_catalog (catalog, 0, i + 20); // for a test, just dump a specific catalog + + setMcal (catalog, i, FALSE, flatcorr); + setMmos (catalog, i, FALSE, flatcorr); + setMgrp (catalog, i, FALSE, flatcorr); + dump_tgroups (catalog, i + 40); + dump_catalog (catalog, 0, i + 40); // for a test, just dump a specific catalog setMgrid (catalog, flatcorr); @@ -147,20 +140,10 @@ plot_chisq (catalog, Ncatalog); } - // if (i < NLOOP - 1) rationalize_mosaics (catalog, Ncatalog); - // if (i % 6 == 1) rationalize_images (); - - // NOTE : in the past, I was not iterating enough before cleaning. make sure we do - // at least 8 loops first -- that should get the systematic errors down to the ~1% - // level, even in cases where we have an even split between photometric data and - // data with 1 mag of extinction. - - if ((i > 8) && (i % 8 == 2)) clean_measures (catalog, Ncatalog, FALSE, flatcorr); - if ((i > 8) && (i % 8 == 3)) clean_stars (catalog, Ncatalog); - if ((i > 8) && (i % 8 == 5)) clean_mosaics (); - if ((i > 8) && (i % 8 == 5)) clean_tgroups (); - if ((i > 8) && (i % 8 == 5)) clean_images (); global_stats (catalog, Ncatalog, flatcorr, i); + + SetZeroPointModes (i); MARKTIME("-- finished loop %d: %f sec\n", i, dtime); + fprintf (stderr, "----- continue loop %d -----\n", i); } @@ -178,7 +161,8 @@ /* set Mcal & Mmos for bad images */ int onlyPoorImages = !CALIBRATE_STACKS_AND_WARPS; - setMcal (catalog, onlyPoorImages, flatcorr); - setMmos (catalog, onlyPoorImages, flatcorr); - setMgrp (catalog, onlyPoorImages, flatcorr); + onlyPoorImages = FALSE; // XXX deactivate for now + setMcal (catalog, 999, onlyPoorImages, flatcorr); // max loop since this is a final pass + setMmos (catalog, 999, onlyPoorImages, flatcorr); // max loop since this is a final pass + setMgrp (catalog, 999, onlyPoorImages, flatcorr); // max loop since this is a final pass MARKTIME("-- finalize Mcal values: %f sec\n", dtime); Index: trunk/Ohana/src/relphot/src/relphot_parallel_images.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_parallel_images.c (revision 41461) +++ trunk/Ohana/src/relphot/src/relphot_parallel_images.c (revision 41462) @@ -114,6 +114,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, i, FALSE, flatcorr); + setMmos (catalog, i, FALSE, 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); @@ -151,6 +151,6 @@ /* 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); + setMcal (catalog, 999, onlyPoorImages, flatcorr); + setMmos (catalog, 999, onlyPoorImages, flatcorr); // setMgrp (catalog, onlyPoorImages, flatcorr); XXX see above MARKTIME("-- finalize Mcal values: %f sec\n", dtime); Index: trunk/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 41461) +++ trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 41462) @@ -9,4 +9,5 @@ # endif +int magStatsByRanking (StatDataSet *dataset, StatType *stats); int magStatsByRankingClipped (StatDataSet *dataset, StatType *stats); int magStatsByRankingIRLS (StatDataSet *dataset, StatType *stats); @@ -17,5 +18,6 @@ int liststats_irls (StatDataSet *dataset, int Npoints, StatType *stats); -# define MAG_STATS_BY_RANKING magStatsByRankingIRLS +// # define MAG_STATS_BY_RANKING magStatsByRankingIRLS +# define MAG_STATS_BY_RANKING magStatsByRanking # define UBERCAL_WEIGHT 100.0 @@ -56,4 +58,6 @@ } +FILE *TESTFILE = NULL; + int setMrel_catalog_alt (Catalog *catalog, int Nc, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) { @@ -67,4 +71,6 @@ SetMrelInfoReset (results); // reset the counters + TESTFILE = fopen ("mytest.dat", "w"); + for (j = 0; j < catalog[Nc].Naverage; j++) { @@ -83,4 +89,6 @@ } } + + fclose (TESTFILE); return (TRUE); } @@ -222,4 +230,5 @@ // ** Choose the calibration (depends on the mode : do I have an image reference or not?) float Mcal = 0.0, Mmos = 0.0, Mgrid = 0.0, Mflat = 0.0, Mgrp = 0.0; + float dMgrp = 0.0; off_t meas = measureOffset + k; if (getImageEntry (meas, cat) < 0) { @@ -235,5 +244,5 @@ Mmos = getMmos (meas, cat); if (isnan(Mmos)) SKIP_THIS_MEAS(Nmos); - Mgrp = getMgrp (meas, cat, measureT[k].airmass); + Mgrp = getMgrp (meas, cat, measureT[k].airmass, &dMgrp); if (isnan(Mgrp)) SKIP_THIS_MEAS(Ngrp); Mgrid = getMgrid (meas, cat); @@ -251,4 +260,7 @@ int useUbercalWeight = isUbercal || (refPhotcode && (code->code == refPhotcode->code)); + // combine the various errors in quadrature here: + float dMsys = hypot (code->photomErrSys, dMgrp); + float Map = NAN; float dMap = NAN; @@ -258,6 +270,6 @@ // NOTE: for PV3, apFluxErr is broken (see pmSourcePhotometry.c:245). we are going to // use PSF flux error instead here: - // dMap = MAX (hypot(measure[k].dMap, code->photomErrSys), MIN_ERROR); - dMap = MAX (hypot(measure[k].dM, code->photomErrSys), MIN_ERROR); + // dMap = MAX (hypot(measure[k].dMap, dMsys), MIN_ERROR); + dMap = MAX (hypot(measure[k].dM, dMsys), MIN_ERROR); if (CHECK_VALID_MAG(Map, dMap)) { @@ -277,7 +289,7 @@ float dMkron; if (isSetMrelFinal) { - dMkron = MAX (hypot(measure[k].dMkron, code->photomErrSys), MIN_ERROR); + dMkron = MAX (hypot(measure[k].dMkron, dMsys), MIN_ERROR); } else { - dMkron = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR); + dMkron = MAX (hypot(measureT[k].dM, dMsys), MIN_ERROR); } if (CHECK_VALID_MAG(Mkron, dMkron)) { @@ -293,5 +305,5 @@ float Mpsf = PhotSysTiny (&measureT[k], &averageT[0], &secfilt[0], MAG_CLASS_PSF); - float dMpsf = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR); + float dMpsf = MAX (hypot(measureT[k].dM, dMsys), MIN_ERROR); if (CHECK_VALID_MAG(Mpsf, dMpsf)) { int Npsf = results->psfData[Nsec].Nlist; @@ -401,4 +413,5 @@ secfilt[Nsec].dMpsfChp = psfstats->error; secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN; + fprintf (TESTFILE, "%d %d : %d : %f %f : %f\n", cat, (int) ave, NrankingPSF, psfstats->mean, psfstats->median, psfstats->mean - psfstats->median); } Index: trunk/Ohana/src/relphot/src/test_fitirls.c =================================================================== --- trunk/Ohana/src/relphot/src/test_fitirls.c (revision 41461) +++ trunk/Ohana/src/relphot/src/test_fitirls.c (revision 41462) @@ -5,8 +5,9 @@ int main (int argc, char **argv) { - int Ntests = 100; - int Npoints = 100; - int Noutliers = 10; - int Nbootstrap = 100; + int Ntests = 100; + int Npoints = 100; + int Noutliers = 10; + int Nbootstrap = 100; + int Niterations = 10; int N; @@ -29,4 +30,9 @@ remove_argument (N, &argc, argv); Nbootstrap = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-Niterations"))) { + remove_argument (N, &argc, argv); + Niterations = atoi (argv[N]); remove_argument (N, &argc, argv); } @@ -88,4 +94,5 @@ FitDataSet dataset; FitDataSetAlloc (&dataset, Npoints + Noutliers, FitSlope, Nbootstrap); + dataset.MaxIterations = Niterations; if (isfinite(ValuePrior) || isfinite(SlopePrior)) { Index: trunk/Ohana/src/relphot/test/relphot.tgroups.dvo =================================================================== --- trunk/Ohana/src/relphot/test/relphot.tgroups.dvo (revision 41461) +++ trunk/Ohana/src/relphot/test/relphot.tgroups.dvo (revision 41462) @@ -451,2 +451,35 @@ set stars_dmag = 0.02 + 0.1 * rnd(tmp) end + +macro plot.test.files + if ($0 != 2) + echo "USAGE: plot.test.files (N)" + break + end + + clear -s + resize 2000 1200 + + section a0 0.00 0.50 0.25 0.50 + section a1 0.25 0.50 0.25 0.50 + section a2 0.50 0.50 0.25 0.50 + section a3 0.75 0.50 0.25 0.50 + + section a4 0.00 0.00 0.25 0.50 + section a5 0.25 0.00 0.25 0.50 + section a6 0.50 0.00 0.25 0.50 + section a7 0.75 0.00 0.25 0.50 + + label -fn helvetica 12 + $BSTYLE = -tickpad 0.2 +xpad 0.2 +ypad 0.2 -ypad 5 + + for i 0 8 + sprintf name tgrp.%05d.%02d.dat $i $1 + data $name + exec ls -l --full-time $name + read Msys 4 Mrel 5 Mcal 7 Mmos 8 Mgrp 9 airmass 11 dMag 12 + + section a$i + lim Msys dMag; box $BSTYLE; plot Msys dMag -pt cir -op 0.2 -c blue + end +end