Index: /branches/eam_branches/relphot.20210521/include/relphot.h =================================================================== --- /branches/eam_branches/relphot.20210521/include/relphot.h (revision 41636) +++ /branches/eam_branches/relphot.20210521/include/relphot.h (revision 41637) @@ -411,5 +411,7 @@ double IMAGE_OFFSET; double NIGHT_SCATTER; -double NIGHT_OFFSET; +double NIGHT_CHISQ; +double MOSAIC_SCATTER; +double MOSAIC_CHISQ; double STAR_SCATTER; double STAR_CHISQ; Index: /branches/eam_branches/relphot.20210521/src/ConfigInit.c =================================================================== --- /branches/eam_branches/relphot.20210521/src/ConfigInit.c (revision 41636) +++ /branches/eam_branches/relphot.20210521/src/ConfigInit.c (revision 41637) @@ -26,11 +26,16 @@ GetConfig (config, "SIGMA_LIM", "%lf", 0, &SIGMA_LIM); + DefConfig ("STAR_SCATTER", "%lf", 0.1, STAR_SCATTER); + DefConfig ("STAR_CHISQ", "%lf", 10.0, STAR_CHISQ); - GetConfig (config, "STAR_SCATTER", "%lf", 0, &STAR_SCATTER); + DefConfig ("NIGHT_SCATTER", "%lf", 0.1, NIGHT_SCATTER); + DefConfig ("NIGHT_CHISQ", "%lf", 10.0, NIGHT_CHISQ); + + DefConfig ("MOSAIC_SCATTER", "%lf", 0.1, MOSAIC_SCATTER); + DefConfig ("MOSAIC_CHISQ", "%lf", 10.0, MOSAIC_CHISQ); GetConfig (config, "IMAGE_SCATTER", "%lf", 0, &IMAGE_SCATTER); GetConfig (config, "IMAGE_OFFSET", "%lf", 0, &IMAGE_OFFSET); - GetConfig (config, "STAR_CHISQ", "%lf", 0, &STAR_CHISQ); GetConfig (config, "GRID_TOOFEW", "%d", 0, &GRID_TOOFEW); GetConfig (config, "STAR_TOOFEW", "%d", 0, &STAR_TOOFEW); @@ -54,6 +59,4 @@ DefConfig ("RELPHOT_IMFIT_SYS_SIGMA_LIM", "%lf", 0.01, IMFIT_SYS_SIGMA_LIM); - DefConfig ("NIGHT_SCATTER", "%lf", 0.05, NIGHT_SCATTER); - DefConfig ("NIGHT_OFFSET", "%lf", 0.50, NIGHT_OFFSET); DefConfig ("SKY_DEPTH", "%d", 2, SKY_DEPTH); Index: /branches/eam_branches/relphot.20210521/src/GridOps.c =================================================================== --- /branches/eam_branches/relphot.20210521/src/GridOps.c (revision 41636) +++ /branches/eam_branches/relphot.20210521/src/GridOps.c (revision 41637) @@ -60,4 +60,6 @@ GridCorr[i] = NULL; } + + fprintf (stderr, "Generating grid corrections for %d photcodes\n", NGridCorr); return; } @@ -102,4 +104,5 @@ result = GridCorr[i]; *Nlast = i; + break; } return result; @@ -127,6 +130,8 @@ Image *images = getimages (&Nimages, NULL); + int NGridReal = 0; + for (int i = 0; i < Nimages; i++) { - int code = images[0].photcode; + int code = images[i].photcode; myAssert (code >= 0, "oops"); myAssert (code < NGridCorr, "oops"); @@ -137,4 +142,6 @@ if (GridCorr[code]) continue; // already created this one + + NGridReal ++; ALLOCATE(GridCorr[code], GridCorrectionType, 1); @@ -157,4 +164,5 @@ } resetMgrid(); // start with values of 0 + fprintf (stderr, "Init grid corrections for %d photcodes\n", NGridReal); } @@ -274,6 +282,5 @@ GridCorr[code]-> Mgrid[ix][iy] = Mgrid; GridCorr[code]->dMgrid[ix][iy] = sqrt(r*(Mgrid2 - Mgrid*Mgrid)); // sample stdev - fprintf (stderr, "grid code %d, %d x %d : %f +/- %f : %d\n", code, ix, iy, - GridCorr[code]-> Mgrid[ix][iy], GridCorr[code]->dMgrid[ix][iy], GridCorr[code]->nMgrid[ix][iy]); + // fprintf (stderr, "grid code %d, %d x %d : %f +/- %f : %d\n", code, ix, iy, GridCorr[code]-> Mgrid[ix][iy], GridCorr[code]->dMgrid[ix][iy], GridCorr[code]->nMgrid[ix][iy]); } } Index: /branches/eam_branches/relphot.20210521/src/MosaicOps.c =================================================================== --- /branches/eam_branches/relphot.20210521/src/MosaicOps.c (revision 41636) +++ /branches/eam_branches/relphot.20210521/src/MosaicOps.c (revision 41637) @@ -1483,6 +1483,8 @@ /* mark mosaic if: abs(Mcal - ) too large, dMcal too large */ -static float MinMaxChiSq = NAN; -static float MinMaxScatter = NAN; +static float MinChiSqLim = NAN; +static float MaxChiSqLim = NAN; +static float MinScatterLim = NAN; +static float MaxScatterLim = NAN; void clean_mosaics () { @@ -1494,4 +1496,7 @@ if (VERBOSE) fprintf (stderr, "marking poor mosaics\n"); + + if (isnan (MaxChiSqLim)) MaxChiSqLim = MOSAIC_CHISQ; + if (isnan (MaxScatterLim)) MaxScatterLim = MOSAIC_SCATTER; ALLOCATE (mlist, double, Nmosaic); @@ -1516,13 +1521,13 @@ liststats (mlist, NULL, NULL, N, &stats); float ChiSqUpper90 = stats.Upper90; - if (isnan (MinMaxChiSq)) MinMaxChiSq = 2.0*stats.median; - float MaxChiSq = MAX (MinMaxChiSq, ChiSqUpper90); + if (isnan (MinChiSqLim)) MinChiSqLim = 2.0*stats.median; // chi-square cut cannot fall below this value (even if this is > MaxChiSqLim) + float ChiSqLimit = MAX(MinChiSqLim, MIN(MaxChiSqLim, ChiSqUpper90)); // chi-square cut should be between MinChiSqLim and MaxChiSqLim liststats (slist, NULL, NULL, N, &stats); float ScatterUpper90 = stats.Upper90; - if (isnan (MinMaxScatter)) MinMaxScatter = 2.0*stats.median; - float MaxScatter = MAX (MinMaxScatter, ScatterUpper90); - - fprintf (stderr, "MOSAIC: Max ChiSq: %f, Max Scatter: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", MaxChiSq, MaxScatter, ChiSqUpper90, ScatterUpper90); + if (isnan (MinScatterLim)) MinScatterLim = 2.0*stats.median; // scatter cut cannto fall below this value (even if this is > MaxScatterLim) + float ScatterLimit = MAX(MinScatterLim, MIN(MaxScatterLim, ScatterUpper90)); + + fprintf (stderr, "MOSAIC: ChiSqLimit: %f, ScatterLimit: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", ChiSqLimit, ScatterLimit, ChiSqUpper90, ScatterUpper90); int Ntotal = 0, Npoor = 0, Nmark = 0, Nscatter = 0, Nfew = 0, Nchisq = 0; @@ -1547,9 +1552,9 @@ Nfew ++; } - if (mosaic[i].stdev > MaxScatter) { + if (mosaic[i].stdev > ScatterLimit) { mark = TRUE; Nscatter ++; } - if (mosaic[i].McalChiSq > MaxChiSq) { + if (mosaic[i].McalChiSq > ChiSqLimit) { mark = TRUE; Nchisq ++; Index: /branches/eam_branches/relphot.20210521/src/StarOps.c =================================================================== --- /branches/eam_branches/relphot.20210521/src/StarOps.c (revision 41636) +++ /branches/eam_branches/relphot.20210521/src/StarOps.c (revision 41637) @@ -576,10 +576,16 @@ } +static float MinChiSqLim = NAN; +static float MaxChiSqLim = NAN; +static float MinScatterLim = NAN; +static float MaxScatterLim = NAN; + void clean_stars (Catalog *catalog, int Ncatalog) { - int i, j, Ndel, Nave, Ntot, mark, Ns, Nscat, Nchi, Nnan; - float dM; - double MaxScatter, MaxChisq; - double *xlist, *slist, *dlist; + int Ndel, Nave, Ntot, Nscat, Nchi, Nnan; + double *xlist, *slist; + + if (isnan (MaxChiSqLim)) MaxChiSqLim = STAR_CHISQ; + if (isnan (MaxScatterLim)) MaxScatterLim = STAR_SCATTER; StatType stats; @@ -589,21 +595,20 @@ /* find Mchisq median -> ChiSq lim must be > median */ - for (i = Ntot = 0; i < Ncatalog; i++) { + for (int i = Ntot = 0; i < Ncatalog; i++) { Ntot += catalog[i].Naverage; } ALLOCATE (xlist, double, Ntot); ALLOCATE (slist, double, Ntot); - ALLOCATE (dlist, double, Ntot); int Nsecfilt = GetPhotcodeNsecfilt (); // eliminate bad stars using the stats for a single secfilt at a time - for (Ns = 0; Ns < Nphotcodes; Ns ++) { + for (int Ns = 0; Ns < Nphotcodes; Ns ++) { int thisCode = photcodes[Ns][0].code; int Nsec = GetPhotcodeNsec(thisCode); - for (i = Ntot = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { + for (int i = Ntot = 0; i < Ncatalog; i++) { + for (int j = 0; j < catalog[i].Naverage; j++) { if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq; @@ -611,32 +616,32 @@ xlist[Ntot] = Mchisq; slist[Ntot] = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp; - dlist[Ntot] = 1; Ntot ++; } } - // XXX the limits for MaxChisq and MaxScatter should be - // user-defined. This uses - liststats (xlist, dlist, NULL, Ntot, &stats); - float ChiSqUpper90 = stats.Upper90; - MaxChisq = MAX (STAR_CHISQ, ChiSqUpper90); - - liststats (slist, dlist, NULL, Ntot, &stats); + liststats (xlist, NULL, NULL, Ntot, &stats); + float ChiSqUpper90 = stats.Upper90; + if (isnan (MinChiSqLim)) MinChiSqLim = 2.0*stats.median; // chi-square cut cannot fall below this value (even if this is > MaxChiSqLim) + float ChiSqLimit = MAX(MinChiSqLim, MIN(MaxChiSqLim, ChiSqUpper90)); // chi-square cut should be between MinChiSqLim and MaxChiSqLim + + liststats (slist, NULL, NULL, Ntot, &stats); float ScatterUpper90 = stats.Upper90; - MaxScatter = MAX (STAR_SCATTER, ScatterUpper90); - fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq); + if (isnan (MinScatterLim)) MinScatterLim = 2.0*stats.median; // scatter cut cannto fall below this value (even if this is > MaxScatterLim) + float ScatterLimit = MAX(MinScatterLim, MIN(MaxScatterLim, ScatterUpper90)); + + fprintf (stderr, "STARS: ChiSqLimit: %f, ScatterLimit: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", ChiSqLimit, ScatterLimit, ChiSqUpper90, ScatterUpper90); Ndel = Nave = Nscat = Nnan = Nchi = 0; - for (i = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { - dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp; + for (int i = 0; i < Ncatalog; i++) { + for (int j = 0; j < catalog[i].Naverage; j++) { + float dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp; float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq; - mark = (dM > MaxScatter) || (isnan(Mchisq)) || (Mchisq > MaxChisq); + int mark = (dM > ScatterLimit) || (isnan(Mchisq)) || (Mchisq > ChiSqLimit); if (mark) { catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_STAR_POOR; Ndel ++; - if (dM > MaxScatter) { Nscat ++; } - if (isnan(Mchisq)) { Nnan ++; } - if (Mchisq > MaxChisq) { Nchi ++; } + if (dM > ScatterLimit) { Nscat ++; } + if (isnan(Mchisq)) { Nnan ++; } + if (Mchisq > ChiSqLimit) { Nchi ++; } } else { catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_SECF_STAR_POOR; @@ -645,9 +650,8 @@ } } - fprintf (stderr, "%d stars marked variable (%d scat, %d nan, %d chi), %d total\n", Ndel, Nscat, Nnan, Nchi, Nave); + fprintf (stderr, "%d of %d stars marked variable (%d scat, %d nan, %d chi)\n", Ndel, Nave, Nscat, Nnan, Nchi); } free (xlist); free (slist); - free (dlist); } Index: /branches/eam_branches/relphot.20210521/src/TGroupOps.c =================================================================== --- /branches/eam_branches/relphot.20210521/src/TGroupOps.c (revision 41636) +++ /branches/eam_branches/relphot.20210521/src/TGroupOps.c (revision 41637) @@ -1265,8 +1265,8 @@ */ -static float MinMaxChiSq = NAN; -static float MaxMaxChiSq = 6.0; -static float MinMaxScatter = NAN; -static float MaxMaxScatter = 0.1; +static float MinChiSqLim = NAN; +static float MaxChiSqLim = NAN; +static float MinScatterLim = NAN; +static float MaxScatterLim = NAN; void clean_tgroups () { @@ -1277,4 +1277,7 @@ if (VERBOSE) fprintf (stderr, "marking poor tgroups\n"); + + if (isnan (MaxChiSqLim)) MaxChiSqLim = NIGHT_CHISQ; + if (isnan (MaxScatterLim)) MaxScatterLim = NIGHT_SCATTER; ALLOCATE (mlist, double, NtgroupTimes*Nphotcodes); @@ -1300,13 +1303,13 @@ liststats (mlist, NULL, NULL, N, &stats); float ChiSqUpper90 = stats.Upper90; - if (isnan (MinMaxChiSq)) MinMaxChiSq = 2.0*stats.median; - float MaxChiSq = MIN(MaxMaxChiSq, MAX (MinMaxChiSq, ChiSqUpper90)); // the upper limit should be between MinMaxChiSq and MaxMaxChiSq + if (isnan (MinChiSqLim)) MinChiSqLim = 2.0*stats.median; // chi-square cut cannot fall below this value (even if this is > MaxChiSqLim) + float ChiSqLimit = MAX(MinChiSqLim, MIN(MaxChiSqLim, ChiSqUpper90)); // chi-square cut should be between MinChiSqLim and MaxChiSqLim liststats (slist, NULL, NULL, N, &stats); float ScatterUpper90 = stats.Upper90; - if (isnan (MinMaxScatter)) MinMaxScatter = 2.0*stats.median; - float MaxScatter = MIN(MaxMaxScatter, MAX (MinMaxScatter, ScatterUpper90)); - - fprintf (stderr, "TGROUPS: Max ChiSq: %f, Max Scatter: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", MaxChiSq, MaxScatter, ChiSqUpper90, ScatterUpper90); + if (isnan (MinScatterLim)) MinScatterLim = 2.0*stats.median; // scatter cut cannto fall below this value (even if this is > MaxScatterLim) + float ScatterLimit = MAX(MinScatterLim, MIN(MaxScatterLim, ScatterUpper90)); + + fprintf (stderr, "TGROUPS: ChiSqLimit: %f, ScatterLimit: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", ChiSqLimit, ScatterLimit, ChiSqUpper90, ScatterUpper90); int Ntotal = 0, Npoor = 0, Nmark = 0, Nscatter = 0, Nchisq = 0, NfewNights = 0, NfewExp = 0; @@ -1337,10 +1340,10 @@ } // scatter too large - if (tgroup[j].stdev > MaxScatter) { + if (tgroup[j].stdev > ScatterLimit) { mark = TRUE; Nscatter ++; } // chisq too large - if (tgroup[j].McalChiSq > MaxChiSq) { + if (tgroup[j].McalChiSq > ChiSqLimit) { mark = TRUE; Nchisq ++; Index: /branches/eam_branches/relphot.20210521/src/initialize.c =================================================================== --- /branches/eam_branches/relphot.20210521/src/initialize.c (revision 41636) +++ /branches/eam_branches/relphot.20210521/src/initialize.c (revision 41637) @@ -60,4 +60,12 @@ fprintf (stderr, "STAR_SCATTER %lf\n", STAR_SCATTER); + fprintf (stderr, "STAR_CHISQ %lf\n", STAR_CHISQ); + + fprintf (stderr, "MOSAIC_SCATTER %lf\n", MOSAIC_SCATTER); + fprintf (stderr, "MOSAIC_CHISQ %lf\n", MOSAIC_CHISQ); + + fprintf (stderr, "NIGHT_SCATTER %lf\n", NIGHT_SCATTER); + fprintf (stderr, "NIGHT_CHISQ %lf\n", NIGHT_CHISQ); + fprintf (stderr, "IMAGE_SCATTER %lf\n", IMAGE_SCATTER); fprintf (stderr, "IMAGE_OFFSET %lf\n", IMAGE_OFFSET); Index: /branches/eam_branches/relphot.20210521/src/relphot_images.c =================================================================== --- /branches/eam_branches/relphot.20210521/src/relphot_images.c (revision 41636) +++ /branches/eam_branches/relphot.20210521/src/relphot_images.c (revision 41637) @@ -55,4 +55,6 @@ if (CALIBRATE_STACKS_AND_WARPS || (NLOOP > 0)) { + initGrid (); // allocate grid correction entries for existing photcodes + /* load catalog data from region files (hostID is 0 since we are not a client */ catalog = load_catalogs (skylist, &Ncatalog, 0, NULL, NULL); @@ -78,6 +80,4 @@ findTGroups (catalog, Ncatalog); MARKTIME("-- set up mosaic indexes: %f sec\n", dtime); - - initGrid (); // allocate grid correction entries for existing photcodes SAVEPLOT = FALSE;