IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 31, 2021, 11:45:18 AM (5 years ago)
Author:
eugene
Message:

major revision of the flagging and resets. cleanup of old cruft and deprecated flag values

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r41624 r41625  
    11# include "relphot.h"
     2
     3# define STAR_BAD (ID_SECF_STAR_FEW | ID_SECF_STAR_POOR)
    24
    35enum {THREAD_RUN, THREAD_DONE};
    … …  
    290292
    291293// setMrel and setMrelOutput are extremely similar, but have slightly different implications:
     294// * both call setMrelCatalog, but setMrel set isSetMrelFinal to FALSE
    292295// * setMrel only uses the internal Tiny structures
    293 // * setMrel skips stars for which there are too few good measurements (setMrelOutput assigns an average mag with only 1)
     296// * setMrel skips stars for which there are too few good measurements
     297// -> (setMrelOutput assigns an average mag for all objects, even those with only 1 measurement)
    294298// * setMrelOutput sets psfQF, psfQFperf, extNsigma, rank info, etc
    295299// * setMrelOutput sets average Map
    296300// * setMrelOutput updates 2MASS average flags
    297301// * setMrelOutput updates average EXT flags (PS1 and 2MASS)
    298 // * setMrelOutput sets average stack & warp photometry (if requested)
     302// * setMrelOutput sets average stack & warp photometry (if requested), setMrel always skips those
     303// * setMrelOutput is NOT threaded it reads/writes catalogs one at a time
    299304
    300305int setMrel (Catalog *catalog, int Ncatalog) {
    … …  
    325330}
    326331
    327 // setMrelOutput is NOT threaded because we read/write catalogs one at a time (for now)
    328332int setMrelOutput (Catalog *catalog, int Ncatalog) {
    329333
    … …  
    443447}
    444448
    445 # define SKIP_THIS_MEAS(REASON) {                               \
    446     catalog[Nc].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM;     \
    447     if (catalog[Nc].measure) {                                  \
    448       catalog[Nc].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;   \
    449     }                                                           \
    450     results->REASON ++;                                         \
    451     continue; }
    452 
    453 
    454449int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {
    455450
    … …  
    471466}
    472467
    473 /* set measure.Mcal for all measures except ID_MEAS_NOCAL and ID_IMAGE_PHOTOM_NOCAL */
     468/* set measure.Mcal for all measures except ID_MEAS_NOCAL */
    474469int setMcalOutput (Catalog *catalog, int Ncatalog) {
    475470
    … …  
    477472  off_t j, k, m;
    478473
    479   MEAS_BAD = ID_MEAS_NOCAL;
    480   IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL;
    481 
    482474  for (i = 0; i < Ncatalog; i++) {
    483475    for (j = 0; j < catalog[i].Naverage; j++) {
    … …  
    485477      m = catalog[i].averageT[j].measureOffset;
    486478      for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
    487         if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
     479        if (catalog[i].measureT[m].dbFlags & ID_MEAS_NOCAL) continue;
    488480        float McalPSF  = getMcal  (m, i, MAG_CLASS_PSF);
    489481        float McalAPER = getMcal  (m, i, MAG_CLASS_APER);
    490482        if (isnan(McalPSF)) continue;
    491         float Mmos  = getMmos  (m, i);
    492         if (isnan(Mmos)) continue;
    493         float Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass, NULL);
    494         if (isnan(Mgrp)) continue;
    495         float Mgrid = getMgridTiny (&catalog[i].measureT[m]);
    496         if (isnan(Mgrid)) continue;
     483
     484        // Mmos, Mgrp have already been transfered to Mcal and should be zero here
     485        // float Mmos  = getMmos  (m, i);
     486        // if (isnan(Mmos)) continue;
     487        // float Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass, NULL);
     488        // if (isnan(Mgrp)) continue;
     489
     490        // Mgrid is applied to Mflat and should be ignored here
     491        // float Mgrid = getMgridTiny (&catalog[i].measureT[m]);
     492        // if (isnan(Mgrid)) continue;
    497493
    498494        // Note that this operation is setting measure->McalAPER to image->McalAPER only
    … …  
    507503
    508504        // set the output calibration
    509         // float Moff = Mmos + Mgrp + Mgrid;
    510         float Moff = Mmos + Mgrid; // Mgrp has already been applied to the image.Mcal
    511         catalog[i].measure[m].McalPSF  = McalPSF + Moff;
    512         catalog[i].measure[m].McalAPER = useStackAper ? McalAPER + Moff : McalPSF + Moff;
     505        catalog[i].measure[m].McalPSF  = McalPSF;
     506        catalog[i].measure[m].McalAPER = useStackAper ? McalAPER : McalPSF;
    513507
    514508        if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
    … …  
    612606    for (i = Ntot = 0; i < Ncatalog; i++) {
    613607      for (j = 0; j < catalog[i].Naverage; j++) {
    614         if ( catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD ) continue;
     608        if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
    615609        float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
    616610        if (isnan(Mchisq)) continue;
    … …  
    621615      }
    622616    }
    623  
     617
     618    // XXX the limits for MaxChisq and MaxScatter should be
     619    // user-defined.  This uses
    624620    liststats (xlist, dlist, NULL, Ntot, &stats);
    625     MaxChisq = MAX (STAR_CHISQ, 2*stats.median);
     621    float ChiSqUpper90 = stats.Upper90;
     622    MaxChisq = MAX (STAR_CHISQ, ChiSqUpper90);
    626623
    627624    liststats (slist, dlist, NULL, Ntot, &stats);
    628     MaxScatter = MAX (STAR_SCATTER, 2*stats.median);
     625    float ScatterUpper90 = stats.Upper90;
     626    MaxScatter = MAX (STAR_SCATTER, ScatterUpper90);
    629627    fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq);
    630628
    … …  
    636634        mark = (dM > MaxScatter) || (isnan(Mchisq)) || (Mchisq > MaxChisq);
    637635        if (mark) {
    638           catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_OBJ_POOR;
     636          catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_STAR_POOR;
    639637          Ndel ++;
    640638          if (dM > MaxScatter)   { Nscat ++; }
    … …  
    642640          if (Mchisq > MaxChisq) { Nchi ++; }
    643641        } else {
    644           catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_OBJ_POOR;
     642          catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_SECF_STAR_POOR;
    645643        }
    646644        Nave ++;
    … …  
    652650  free (slist);
    653651  free (dlist);
    654 }
    655 
    656 // clean_measures examines the stats for a single star.  It measures the INNER 50% mean
    657 // and sigma, it then re-measures the mean and sigma using all stars with NSIGMA_CLIP (3)
    658 // sigma of the INNER 50% mean.  it then flags any measurements which are more than
    659 // NSIGMA_REJECT (5) sigma of the mean
    660 
    661 // this function only operations on the PSF magnitudes
    662 
    663 // XXX EAM clean_measures is no longer used -- we rely on IRLS to down-weight outliers
    664 
    665 # define NSIGMA_CLIP 3.0
    666 # define NSIGMA_REJECT 5.0
    667 void clean_measures (Catalog *catalog, int Ncatalog, int final) {
    668 
    669   off_t j, k, m, Nmax, Ndel, Nave;
    670   int i, N, image_bad, TOOFEW;
    671   off_t *ilist;
    672   double *tlist, *list, *dlist;
    673   float Msys, Mcal, Mmos, Mgrp, Mgrid;
    674   int Ncal, Nmos, Ngrp, Ngrid, Nfew;
    675 
    676   int Nsecfilt = GetPhotcodeNsecfilt ();
    677 
    678   if (VERBOSE) fprintf (stderr, "marking poor measures\n");
    679   /* Nmeasure is now different, need to reallocate */
    680   Nmax = 0;
    681   for (i = 0; i < Ncatalog; i++) {
    682     for (j = 0; j < catalog[i].Naverage; j++) {
    683       Nmax = MAX (Nmax, catalog[i].averageT[j].Nmeasure);
    684     }
    685   }
    686   ALLOCATE (list, double, MAX (1, Nmax));
    687   ALLOCATE (dlist, double, MAX (1, Nmax));
    688   ALLOCATE (ilist, off_t, Nmax);
    689   ALLOCATE (tlist, double, Nmax);
    690  
    691   /* it makes no sense to mark 3-sigma outliers with <5 measurements */
    692   TOOFEW = MAX (5, STAR_TOOFEW);
    693 
    694   // stats structures for inner and full stats
    695   StatType instats, stats;
    696   liststats_setmode (&instats, "INNER_MEAN");
    697   liststats_setmode (&stats, "MEAN");
    698 
    699   Ndel = Nave = 0;
    700   Ncal = Nmos = Ngrp = Ngrid = Nfew = 0;
    701   for (i = 0; i < Ncatalog; i++) {
    702     for (j = 0; j < catalog[i].Naverage; j++) {
    703 
    704       int Ns;
    705       for (Ns = 0; Ns < Nphotcodes; Ns++) {
    706        
    707         int thisCode = photcodes[Ns][0].code;
    708         int Nsec = GetPhotcodeNsec(thisCode);
    709        
    710         /* on final processing, skip stars already measured */
    711         if (final && catalog[i].found_t[Nsecfilt*j + Nsec]) continue; 
    712 
    713         /* skip bad stars to prevent them from becoming good (on inner sample) */
    714         if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; 
    715 
    716         /* accumulate list of valid measurements */
    717         m = catalog[i].averageT[j].measureOffset;
    718         N = 0;
    719         for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
    720           // skip measurements that do not match the current photcode
    721           int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
    722           if (ecode != thisCode) { continue; }
    723 
    724           // SKIP gpc1 stack data
    725           if (isGPC1stack(catalog[i].measureT[m].photcode)) continue;
    726          
    727           // SKIP gpc1 forced-warp data
    728           if (isGPC1warp(catalog[i].measureT[m].photcode)) continue;
    729          
    730           // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
    731           Mcal  = getMcal  (m, i, MAG_CLASS_PSF);
    732           if (isnan(Mcal)) { Ncal ++; continue; }
    733           Mmos  = getMmos  (m, i);
    734           if (isnan(Mmos)) { Nmos ++; continue; }
    735           Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass, NULL);
    736           if (isnan(Mgrp)) { Ngrp ++; continue; }
    737           Mgrid = getMgridTiny (&catalog[i].measureT[m]);
    738           if (isnan(Mgrid)) { Ngrid ++; continue; }
    739 
    740           Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF);
    741           list[N] = Msys - Mcal - Mmos - Mgrp - Mgrid;
    742           dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR);
    743           N++;
    744         }
    745         if (N <= TOOFEW) {
    746           Nfew ++;
    747           continue;
    748         }
    749 
    750         /* 3-sigma clip based on stats of inner 50% */
    751 
    752         // calculated mean of inner 50%
    753         liststats (list, dlist, NULL, N, &instats);
    754         instats.sigma = MAX (MIN_ERROR, instats.sigma); /* if measurements agree too well, sigma -> 0.0 */
    755 
    756         // ignore entries > 3sigma from inner mean
    757         for (k = m = 0; k < N; k++) {
    758           if (fabs (list[k] - instats.median) < NSIGMA_CLIP*instats.sigma) {
    759             list[m] = list[k];
    760             m++;
    761           }
    762         }
    763         // recalculate the mean & sigma of the accepted measurements
    764         liststats (list, dlist, NULL, m, &stats);
    765         stats.sigma = MAX (MIN_ERROR, stats.sigma);
    766 
    767         /* apply to list of all relevant measurements, including IMAGE_POOR & IMAGE_FEW */
    768         image_bad = IMAGE_BAD;
    769         IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL;
    770         m = catalog[i].averageT[j].measureOffset;
    771         N = 0;
    772         for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
    773           // skip measurements that do not match the current photcode
    774           int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
    775           if (ecode != thisCode) { continue; }
    776 
    777           // SKIP gpc1 stack data
    778           if (isGPC1stack(catalog[i].measureT[m].photcode)) continue;
    779          
    780           // SKIP gpc1 forced-warp data
    781           if (isGPC1warp(catalog[i].measureT[m].photcode)) continue;
    782          
    783           // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
    784           Mcal  = getMcal  (m, i, MAG_CLASS_PSF);
    785           if (isnan(Mcal)) continue;
    786           Mmos  = getMmos  (m, i);
    787           if (isnan(Mmos)) continue;
    788           Mgrid = getMgridTiny (&catalog[i].measureT[m]);
    789           if (isnan(Mgrid)) continue;
    790 
    791           Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF);
    792           list[N] = Msys - Mcal - Mmos - Mgrid;
    793           dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR);
    794           ilist[N] = m;
    795           N++;
    796           Nave ++;
    797         }
    798         if (N < TOOFEW) continue;
    799 
    800         /* mark bad measures (> 5 sigma deviant) */
    801         for (k = 0; k < N; k++) {
    802           // treat the scatter of the star as a systematic term.  this is a bit of an
    803           // over-estimage (a perfect Gauss distribution with perfect errors would have
    804           // mySigma = sqrt(2) too large)
    805           float mySigma = hypot (stats.sigma, dlist[k]);
    806           if (fabs (list[k] - stats.median) > NSIGMA_REJECT*mySigma) {
    807             catalog[i].measureT[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM;
    808             if (final) {
    809               // for the final pass, we have a duplicate set of values in measure and measureT
    810               catalog[i].measure[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM;
    811             }
    812             Ndel ++;
    813           } else {
    814             catalog[i].measureT[ilist[k]].dbFlags &= ~ID_MEAS_POOR_PHOTOM;
    815             if (final) {
    816               // for the final pass, we have a duplicate set of values in measure and measureT
    817               catalog[i].measure[ilist[k]].dbFlags &= ~ID_MEAS_POOR_PHOTOM;
    818             }
    819           }
    820         }
    821         IMAGE_BAD = image_bad;
    822       }
    823     }
    824   }
    825   if (VERBOSE) fprintf (stderr, OFF_T_FMT" measures marked poor, "OFF_T_FMT" total\n", Ndel, Nave);
    826 
    827   free (list);
    828   free (dlist);
    829   free (ilist);
    830   free (tlist);
    831652}
    832653
Note: See TracChangeset for help on using the changeset viewer.