IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 41637


Ignore:
Timestamp:
Jun 4, 2021, 10:46:36 AM (5 years ago)
Author:
eugene
Message:

add user-control over night, mosaic, image chisq & scatter cuts

Location:
branches/eam_branches/relphot.20210521
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/relphot.20210521/include/relphot.h

    r41625 r41637  
    411411double IMAGE_OFFSET;
    412412double NIGHT_SCATTER;
    413 double NIGHT_OFFSET;
     413double NIGHT_CHISQ;
     414double MOSAIC_SCATTER;
     415double MOSAIC_CHISQ;
    414416double STAR_SCATTER;
    415417double STAR_CHISQ;

  • branches/eam_branches/relphot.20210521/src/ConfigInit.c

    r41603 r41637  
    2626  GetConfig (config, "SIGMA_LIM",              "%lf", 0, &SIGMA_LIM);
    2727
     28  DefConfig ("STAR_SCATTER",                    "%lf",  0.1, STAR_SCATTER);
     29  DefConfig ("STAR_CHISQ",                      "%lf", 10.0, STAR_CHISQ);
    2830
    29   GetConfig (config, "STAR_SCATTER",           "%lf", 0, &STAR_SCATTER);
     31  DefConfig ("NIGHT_SCATTER",                   "%lf",  0.1, NIGHT_SCATTER);
     32  DefConfig ("NIGHT_CHISQ",                     "%lf", 10.0, NIGHT_CHISQ);
     33
     34  DefConfig ("MOSAIC_SCATTER",                  "%lf",  0.1, MOSAIC_SCATTER);
     35  DefConfig ("MOSAIC_CHISQ",                    "%lf", 10.0, MOSAIC_CHISQ);
    3036
    3137  GetConfig (config, "IMAGE_SCATTER",          "%lf", 0, &IMAGE_SCATTER);
    3238  GetConfig (config, "IMAGE_OFFSET",           "%lf", 0, &IMAGE_OFFSET);
    3339
    34   GetConfig (config, "STAR_CHISQ",             "%lf", 0, &STAR_CHISQ);
    3540  GetConfig (config, "GRID_TOOFEW",            "%d",  0, &GRID_TOOFEW);
    3641  GetConfig (config, "STAR_TOOFEW",            "%d",  0, &STAR_TOOFEW);
    … …  
    5459
    5560  DefConfig ("RELPHOT_IMFIT_SYS_SIGMA_LIM", "%lf", 0.01, IMFIT_SYS_SIGMA_LIM);
    56   DefConfig ("NIGHT_SCATTER",               "%lf", 0.05, NIGHT_SCATTER);
    57   DefConfig ("NIGHT_OFFSET",                "%lf", 0.50, NIGHT_OFFSET);
    5861  DefConfig ("SKY_DEPTH",                   "%d",     2, SKY_DEPTH);
    5962

  • branches/eam_branches/relphot.20210521/src/GridOps.c

    r41632 r41637  
    6060    GridCorr[i] = NULL;
    6161  }
     62
     63  fprintf (stderr, "Generating grid corrections for %d photcodes\n", NGridCorr);
    6264  return;
    6365}
    … …  
    102104    result = GridCorr[i];
    103105    *Nlast = i;
     106    break;
    104107  }
    105108  return result;
    … …  
    127130  Image *images = getimages (&Nimages, NULL);
    128131
     132  int NGridReal = 0;
     133
    129134  for (int i = 0; i < Nimages; i++) {
    130     int code = images[0].photcode;
     135    int code = images[i].photcode;
    131136    myAssert (code >= 0, "oops");
    132137    myAssert (code < NGridCorr, "oops");
    … …  
    137142
    138143    if (GridCorr[code]) continue; // already created this one
     144
     145    NGridReal ++;
    139146
    140147    ALLOCATE(GridCorr[code], GridCorrectionType, 1);
    … …  
    157164  }
    158165  resetMgrid(); // start with values of 0
     166  fprintf (stderr, "Init grid corrections for %d photcodes\n", NGridReal);
    159167}
    160168
    … …  
    274282        GridCorr[code]-> Mgrid[ix][iy] = Mgrid;
    275283        GridCorr[code]->dMgrid[ix][iy] = sqrt(r*(Mgrid2 - Mgrid*Mgrid)); // sample stdev
    276         fprintf (stderr, "grid code %d, %d x %d : %f +/- %f : %d\n", code, ix, iy,
    277                  GridCorr[code]-> Mgrid[ix][iy], GridCorr[code]->dMgrid[ix][iy], GridCorr[code]->nMgrid[ix][iy]);
     284        // 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]);
    278285      }
    279286    }

  • branches/eam_branches/relphot.20210521/src/MosaicOps.c

    r41625 r41637  
    14831483/* mark mosaic if: abs(Mcal - <Mcal>) too large, dMcal too large */
    14841484
    1485 static float MinMaxChiSq = NAN;
    1486 static float MinMaxScatter = NAN;
     1485static float MinChiSqLim = NAN;
     1486static float MaxChiSqLim = NAN;
     1487static float MinScatterLim = NAN;
     1488static float MaxScatterLim = NAN;
    14871489
    14881490void clean_mosaics () {
    … …  
    14941496
    14951497  if (VERBOSE) fprintf (stderr, "marking poor mosaics\n");
     1498
     1499  if (isnan (MaxChiSqLim))   MaxChiSqLim   = MOSAIC_CHISQ;
     1500  if (isnan (MaxScatterLim)) MaxScatterLim = MOSAIC_SCATTER;
    14961501
    14971502  ALLOCATE (mlist, double, Nmosaic);
    … …  
    15161521  liststats (mlist, NULL, NULL, N, &stats);
    15171522  float ChiSqUpper90 = stats.Upper90;
    1518   if (isnan (MinMaxChiSq)) MinMaxChiSq = 2.0*stats.median;
    1519   float MaxChiSq = MAX (MinMaxChiSq, ChiSqUpper90);
     1523  if (isnan (MinChiSqLim)) MinChiSqLim = 2.0*stats.median;             // chi-square cut cannot fall below this value (even if this is > MaxChiSqLim)
     1524  float ChiSqLimit = MAX(MinChiSqLim, MIN(MaxChiSqLim, ChiSqUpper90)); // chi-square cut should be between MinChiSqLim and MaxChiSqLim
    15201525
    15211526  liststats (slist, NULL, NULL, N, &stats);
    15221527  float ScatterUpper90 = stats.Upper90;
    1523   if (isnan (MinMaxScatter)) MinMaxScatter = 2.0*stats.median;
    1524   float MaxScatter = MAX (MinMaxScatter, ScatterUpper90);
    1525 
    1526   fprintf (stderr, "MOSAIC: Max ChiSq: %f, Max Scatter: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", MaxChiSq, MaxScatter, ChiSqUpper90, ScatterUpper90);
     1528  if (isnan (MinScatterLim)) MinScatterLim = 2.0*stats.median;         // scatter cut cannto fall below this value (even if this is > MaxScatterLim)
     1529  float ScatterLimit = MAX(MinScatterLim, MIN(MaxScatterLim, ScatterUpper90));
     1530
     1531  fprintf (stderr, "MOSAIC: ChiSqLimit: %f, ScatterLimit: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", ChiSqLimit, ScatterLimit, ChiSqUpper90, ScatterUpper90);
    15271532 
    15281533  int Ntotal = 0, Npoor = 0, Nmark = 0, Nscatter = 0, Nfew = 0, Nchisq = 0;
    … …  
    15471552      Nfew ++;
    15481553    }
    1549     if (mosaic[i].stdev > MaxScatter) {
     1554    if (mosaic[i].stdev > ScatterLimit) {
    15501555      mark = TRUE;
    15511556      Nscatter ++;
    15521557    }
    1553     if (mosaic[i].McalChiSq > MaxChiSq) {
     1558    if (mosaic[i].McalChiSq > ChiSqLimit) {
    15541559      mark = TRUE;
    15551560      Nchisq ++;

  • branches/eam_branches/relphot.20210521/src/StarOps.c

    r41625 r41637  
    576576}
    577577
     578static float MinChiSqLim = NAN;
     579static float MaxChiSqLim = NAN;
     580static float MinScatterLim = NAN;
     581static float MaxScatterLim = NAN;
     582
    578583void clean_stars (Catalog *catalog, int Ncatalog) {
    579584
    580   int i, j, Ndel, Nave, Ntot, mark, Ns, Nscat, Nchi, Nnan;
    581   float dM;
    582   double MaxScatter, MaxChisq;
    583   double *xlist, *slist, *dlist;
     585  int Ndel, Nave, Ntot, Nscat, Nchi, Nnan;
     586  double *xlist, *slist;
     587
     588  if (isnan (MaxChiSqLim))   MaxChiSqLim   = STAR_CHISQ;
     589  if (isnan (MaxScatterLim)) MaxScatterLim = STAR_SCATTER;
    584590
    585591  StatType stats;
    … …  
    589595
    590596  /* find Mchisq median -> ChiSq lim must be > median */
    591   for (i = Ntot = 0; i < Ncatalog; i++) {
     597  for (int i = Ntot = 0; i < Ncatalog; i++) {
    592598    Ntot += catalog[i].Naverage;
    593599  }
    594600  ALLOCATE (xlist, double, Ntot);
    595601  ALLOCATE (slist, double, Ntot);
    596   ALLOCATE (dlist, double, Ntot);
    597602
    598603  int Nsecfilt = GetPhotcodeNsecfilt ();
    599604
    600605  // eliminate bad stars using the stats for a single secfilt at a time
    601   for (Ns = 0; Ns < Nphotcodes; Ns ++) {
     606  for (int Ns = 0; Ns < Nphotcodes; Ns ++) {
    602607   
    603608    int thisCode = photcodes[Ns][0].code;
    604609    int Nsec = GetPhotcodeNsec(thisCode);
    605610
    606     for (i = Ntot = 0; i < Ncatalog; i++) {
    607       for (j = 0; j < catalog[i].Naverage; j++) {
     611    for (int i = Ntot = 0; i < Ncatalog; i++) {
     612      for (int j = 0; j < catalog[i].Naverage; j++) {
    608613        if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
    609614        float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
    … …  
    611616        xlist[Ntot] = Mchisq;
    612617        slist[Ntot] = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp;
    613         dlist[Ntot] = 1;
    614618        Ntot ++;
    615619      }
    616620    }
    617621
    618     // XXX the limits for MaxChisq and MaxScatter should be
    619     // user-defined.  This uses
    620     liststats (xlist, dlist, NULL, Ntot, &stats);
    621     float ChiSqUpper90 = stats.Upper90;
    622     MaxChisq = MAX (STAR_CHISQ, ChiSqUpper90);
    623 
    624     liststats (slist, dlist, NULL, Ntot, &stats);
     622    liststats (xlist, NULL, NULL, Ntot, &stats);
     623    float ChiSqUpper90 = stats.Upper90;
     624    if (isnan (MinChiSqLim)) MinChiSqLim = 2.0*stats.median;             // chi-square cut cannot fall below this value (even if this is > MaxChiSqLim)
     625    float ChiSqLimit = MAX(MinChiSqLim, MIN(MaxChiSqLim, ChiSqUpper90)); // chi-square cut should be between MinChiSqLim and MaxChiSqLim
     626
     627    liststats (slist, NULL, NULL, Ntot, &stats);
    625628    float ScatterUpper90 = stats.Upper90;
    626     MaxScatter = MAX (STAR_SCATTER, ScatterUpper90);
    627     fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq);
     629    if (isnan (MinScatterLim)) MinScatterLim = 2.0*stats.median;         // scatter cut cannto fall below this value (even if this is > MaxScatterLim)
     630    float ScatterLimit = MAX(MinScatterLim, MIN(MaxScatterLim, ScatterUpper90));
     631
     632    fprintf (stderr, "STARS: ChiSqLimit: %f, ScatterLimit: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", ChiSqLimit, ScatterLimit, ChiSqUpper90, ScatterUpper90);
    628633
    629634    Ndel = Nave = Nscat = Nnan = Nchi = 0;
    630     for (i = 0; i < Ncatalog; i++) {
    631       for (j = 0; j < catalog[i].Naverage; j++) {
    632         dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp;
     635    for (int i = 0; i < Ncatalog; i++) {
     636      for (int j = 0; j < catalog[i].Naverage; j++) {
     637        float dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp;
    633638        float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
    634         mark = (dM > MaxScatter) || (isnan(Mchisq)) || (Mchisq > MaxChisq);
     639        int mark = (dM > ScatterLimit) || (isnan(Mchisq)) || (Mchisq > ChiSqLimit);
    635640        if (mark) {
    636641          catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_STAR_POOR;
    637642          Ndel ++;
    638           if (dM > MaxScatter)   { Nscat ++; }
    639           if (isnan(Mchisq))     { Nnan ++; }
    640           if (Mchisq > MaxChisq) { Nchi ++; }
     643          if (dM > ScatterLimit)   { Nscat ++; }
     644          if (isnan(Mchisq))       { Nnan ++; }
     645          if (Mchisq > ChiSqLimit) { Nchi ++; }
    641646        } else {
    642647          catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_SECF_STAR_POOR;
    … …  
    645650      }
    646651    }
    647     fprintf (stderr, "%d stars marked variable (%d scat, %d nan, %d chi), %d total\n", Ndel, Nscat, Nnan, Nchi, Nave);
     652    fprintf (stderr, "%d of %d stars marked variable (%d scat, %d nan, %d chi)\n", Ndel, Nave, Nscat, Nnan, Nchi);
    648653  }
    649654  free (xlist);
    650655  free (slist);
    651   free (dlist);
    652656}
    653657

  • branches/eam_branches/relphot.20210521/src/TGroupOps.c

    r41625 r41637  
    12651265 */
    12661266
    1267 static float MinMaxChiSq = NAN;
    1268 static float MaxMaxChiSq = 6.0;
    1269 static float MinMaxScatter = NAN;
    1270 static float MaxMaxScatter = 0.1;
     1267static float MinChiSqLim = NAN;
     1268static float MaxChiSqLim = NAN;
     1269static float MinScatterLim = NAN;
     1270static float MaxScatterLim = NAN;
    12711271
    12721272void clean_tgroups () {
    … …  
    12771277
    12781278  if (VERBOSE) fprintf (stderr, "marking poor tgroups\n");
     1279
     1280  if (isnan (MaxChiSqLim))   MaxChiSqLim   = NIGHT_CHISQ;
     1281  if (isnan (MaxScatterLim)) MaxScatterLim = NIGHT_SCATTER;
    12791282
    12801283  ALLOCATE (mlist, double, NtgroupTimes*Nphotcodes);
    … …  
    13001303  liststats (mlist, NULL, NULL, N, &stats);
    13011304  float ChiSqUpper90 = stats.Upper90;
    1302   if (isnan (MinMaxChiSq)) MinMaxChiSq = 2.0*stats.median;
    1303   float MaxChiSq = MIN(MaxMaxChiSq, MAX (MinMaxChiSq, ChiSqUpper90)); // the upper limit should be between MinMaxChiSq and MaxMaxChiSq
     1305  if (isnan (MinChiSqLim)) MinChiSqLim = 2.0*stats.median;             // chi-square cut cannot fall below this value (even if this is > MaxChiSqLim)
     1306  float ChiSqLimit = MAX(MinChiSqLim, MIN(MaxChiSqLim, ChiSqUpper90)); // chi-square cut should be between MinChiSqLim and MaxChiSqLim
    13041307
    13051308  liststats (slist, NULL, NULL, N, &stats);
    13061309  float ScatterUpper90 = stats.Upper90;
    1307   if (isnan (MinMaxScatter)) MinMaxScatter = 2.0*stats.median;
    1308   float MaxScatter = MIN(MaxMaxScatter, MAX (MinMaxScatter, ScatterUpper90));
    1309 
    1310   fprintf (stderr, "TGROUPS: Max ChiSq: %f, Max Scatter: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", MaxChiSq, MaxScatter, ChiSqUpper90, ScatterUpper90);
     1310  if (isnan (MinScatterLim)) MinScatterLim = 2.0*stats.median;         // scatter cut cannto fall below this value (even if this is > MaxScatterLim)
     1311  float ScatterLimit = MAX(MinScatterLim, MIN(MaxScatterLim, ScatterUpper90));
     1312
     1313  fprintf (stderr, "TGROUPS: ChiSqLimit: %f, ScatterLimit: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", ChiSqLimit, ScatterLimit, ChiSqUpper90, ScatterUpper90);
    13111314 
    13121315  int Ntotal = 0, Npoor = 0, Nmark = 0, Nscatter = 0, Nchisq = 0, NfewNights = 0, NfewExp = 0;
    … …  
    13371340      }
    13381341      // scatter too large
    1339       if (tgroup[j].stdev > MaxScatter) {
     1342      if (tgroup[j].stdev > ScatterLimit) {
    13401343        mark = TRUE;
    13411344        Nscatter ++;
    13421345      }
    13431346      // chisq too large
    1344       if (tgroup[j].McalChiSq > MaxChiSq) {
     1347      if (tgroup[j].McalChiSq > ChiSqLimit) {
    13451348        mark = TRUE;
    13461349        Nchisq ++;

  • branches/eam_branches/relphot.20210521/src/initialize.c

    r41625 r41637  
    6060
    6161    fprintf (stderr, "STAR_SCATTER           %lf\n", STAR_SCATTER);
     62    fprintf (stderr, "STAR_CHISQ             %lf\n", STAR_CHISQ);
     63
     64    fprintf (stderr, "MOSAIC_SCATTER         %lf\n", MOSAIC_SCATTER);
     65    fprintf (stderr, "MOSAIC_CHISQ           %lf\n", MOSAIC_CHISQ);
     66
     67    fprintf (stderr, "NIGHT_SCATTER          %lf\n", NIGHT_SCATTER);
     68    fprintf (stderr, "NIGHT_CHISQ            %lf\n", NIGHT_CHISQ);
     69
    6270    fprintf (stderr, "IMAGE_SCATTER          %lf\n", IMAGE_SCATTER);
    6371    fprintf (stderr, "IMAGE_OFFSET           %lf\n", IMAGE_OFFSET);

  • branches/eam_branches/relphot.20210521/src/relphot_images.c

    r41631 r41637  
    5555
    5656  if (CALIBRATE_STACKS_AND_WARPS || (NLOOP > 0)) {
     57    initGrid (); // allocate grid correction entries for existing photcodes
     58
    5759    /* load catalog data from region files (hostID is 0 since we are not a client */
    5860    catalog = load_catalogs (skylist, &Ncatalog, 0, NULL, NULL);
    … …  
    7880    findTGroups (catalog, Ncatalog);
    7981    MARKTIME("-- set up mosaic indexes: %f sec\n", dtime);
    80 
    81     initGrid (); // allocate grid correction entries for existing photcodes
    8282
    8383    SAVEPLOT = FALSE;
Note: See TracChangeset for help on using the changeset viewer.