IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38466 for trunk


Ignore:
Timestamp:
Jun 12, 2015, 10:28:36 PM (11 years ago)
Author:
eugene
Message:

merge changes from EAM branch ipp-20150419

Location:
trunk/Ohana/src/relphot
Files:
16 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relphot/include/relphot.h

    r38441 r38466  
    9696
    9797typedef struct {
     98  double *flxlist;           // list of measure.mag values for a given star
     99  double *errlist;            // mag errors for a star
     100  double *wgtlist;            // weights to use for mean mags
     101  int    *ranking;            // weights to use for mean mags
     102  int    *measSeq;            // weights to use for mean mags
     103  int     Nlist;
     104} StatDataSet;
     105
     106typedef struct {
    98107  int Nfew;
    99108  int Ncode;
    … …  
    105114
    106115  // NOTE: the following arrays are (possibly) pre-allocated and carried down to each
    107   // thread.  The first 3 (list, dlist, wlist) are used in all relphot analyses; the
    108   // others are only used on the final output steps.
    109 
    110   double *Mpsflist;           // list of measure.mag values for a given star
    111   double *dpsflist;           // mag errors for a star
    112   double *wpsflist;           // weights to use for mean mags
    113 
    114   double *Maplist;            // ap mags for a star
    115   double *daplist;            // ap mags for a star
    116   double *waplist;            // ap mags for a star
    117 
    118   double *Mkronlist;          // kron mags for a star
    119   double *dkronlist;          // kron mag errors
    120   double *wkronlist;          // kron mag errors
     116  // thread.  The psfData are used in all relphot analyses; the others are only used on
     117  // the final output steps.
     118
     119  int Nsecfilt;
     120
     121  StatDataSet  *psfData; // one is allocated for each primary (average) photcode
     122  StatDataSet *aperData;
     123  StatDataSet *kronData;
    121124
    122125  double *psfqf_list;         // psfqf for all filters
    123126  double *psfqfperf_list;     // psfqfperf for all filters
    124127  double *stargal_list;       // stargal for all filters
     128
     129  int   *havePS1;          // this secfilt has synthetic mags
     130  int   *haveSYN;          // this secfilt has synthetic mags
     131  int   *measSYN;          // this measurement is the synthetic mag for this secfilt
     132  int   *needSYN;          // this secfilt mag should use synthetic mags
     133  float *minSYN;           // minimum synthetic mag below which synthetic should be forced
     134 
     135  float *psfQfMax;            // max psfQf value for this secfilt
     136  float *psfQfPerfMax;        // max psfQfperf value for this secfilt
     137
     138  int *Nmeas;                 // count of PS1 exposure (chip) measurements for this secfilt
     139  int *NmeasGood;             // count of PS1 exposure (chip) measurements for this secfilt
     140  int *Next;                  // count of PS1 exposure (chip) measurements for this secfilt
     141  int *NexpPS1;               // count of PS1 exposure (chip) measurements for this secfilt
     142  int *haveUbercal;           // does this secfilt have any ubercal data?
     143
     144  int *tessID;                // tess,proj,skycell to use for warp and diff analysis
     145  int *projID;
     146  int *skycellID;
     147
     148  float *minUbercalDist;
    125149
    126150  StatType psfstats;
    … …  
    375399void          liststats_setmode   PROTO((StatType *stats, char *strmode));
    376400int           liststats           PROTO((double *value, double *dvalue, double *wvalue, int N, StatType *stats));
     401int           liststats_init      PROTO((StatType *stats));
    377402Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath, char *syncfile));
    378403Catalog      *load_catalogs_parallel PROTO((SkyList *sky, int *Ncatalog, char *syncfile));
    … …  
    412437int           setMrel             PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
    413438void          setMrelFinal        PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage));
    414 int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr));
     439int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
    415440int           setMave             PROTO((Catalog *catalog, int Ncatalog));
    416441void          set_ZP              PROTO((double ZERO));
    … …  
    474499
    475500int print_measure_set_alt (Average *average, SecFilt *secfilt, Measure *measure);
    476 int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt);
    477 int setMrelAverageExposure (off_t meas, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found);
    478 int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t Nmeasure);
    479 int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found, off_t Nmeasure);
     501int setMrel_catalog_alt (Catalog *catalog, int Nc, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt);
     502int setMrelAverageExposure (Catalog *catalog, int cat, off_t ave, int Nsecfilt, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results);
     503int setMrelAverageStack (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr);
     504int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results);
    480505
    481506int setGlobalObjStats (Average *average, Measure *measure);
    482507
    483508void SetMrelInfoInit (SetMrelInfo *results, int allocLists);
     509void SetMrelInfoFree (SetMrelInfo *results);
     510void SetMrelInfoReset (SetMrelInfo *results);
     511void SetMrelInfoResetObject (SetMrelInfo *results);
     512
     513void StatDataSetFree (StatDataSet *dataset, int Nsecfilt);
     514StatDataSet *StatDataSetAlloc (int Nsecfilt, int Nmax);
    484515
    485516int init_synthetic_mags ();
    … …  
    533564int isGPC1synth (int photcode);
    534565int whichGPC1filter (int photcode);
    535 
    536 
    537 SynthZeroPoints *SynthZeroPointsLoad (char *filename);
     566int is2MASS (int photcode);
     567int isTYCHO (int photcode);
     568
     569int magStatsByRanking (StatDataSet *dataset, StatType *stats);
     570
     571int SynthZeroPointsLoad (char *filename);
     572SynthZeroPoints *SynthZeroPointsGet ();
    538573
    539574int relphot_synthphot (int hostID, char *hostpath);
    … …  
    542577int relphot_synthphot_catalog (Catalog *catalog, SynthZeroPoints *zpts);
    543578int relphot_synthphot_average (Average *average, SecFilt *secfilt, Measure *measure, SynthZeroPoints *zpts);
     579
     580void setMeasureRank (Catalog *catalog);
     581int getImageFlags (off_t meas, int cat);
     582int getMosaicFlags (off_t meas, int cat);
     583

  • trunk/Ohana/src/relphot/src/ImageOps.c

    r37807 r38466  
    397397  i = MeasureToImage[cat][meas];
    398398  return (i);
     399}
     400
     401int getImageFlags (off_t meas, int cat) {
     402
     403  off_t i;
     404
     405  if (!MeasureToImage) return 0;
     406
     407  i = MeasureToImage[cat][meas];
     408
     409  return (image[i].flags);
    399410}
    400411
    … …  
    681692      list[N] = Msys - Mrel - Mmos - Mgrid + Mflat;
    682693      dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
     694      if (fabs(list[N]) > 0.03) {
     695        // fprintf (stderr, "deviant\n");
     696      }
     697
    683698      if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
    684699        Mlist[Nbright] = list[N];
    … …  
    711726
    712727    if (PLOTSTUFF) {
    713       fprintf (stderr, "Mcal for : %s\n", image[i].name);
     728      fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].Mcal, image[i].dMcal);
    714729      plot_setMcal (list, N, &stats, CLOUD_TOLERANCE);
    715730    }

  • trunk/Ohana/src/relphot/src/MosaicOps.c

    r37807 r38466  
    882882}
    883883
     884int getMosaicFlags (off_t meas, int cat) {
     885
     886  if (!MOSAIC_ZEROPT) return (0);
     887
     888  // unassigned measurements belong to simple chips
     889  int i = MeasureToMosaic[cat][meas];
     890  if (i == -1) return (0);
     891
     892  return (mosaic[i].flags);
     893}
     894
    884895typedef struct {
    885896  int Nfew;

  • trunk/Ohana/src/relphot/src/StarOps.c

    r38441 r38466  
    11# include "relphot.h"
    2 
    3 static int Nmax;
    42
    53enum {THREAD_RUN, THREAD_DONE};
    … …  
    2018int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure);
    2119
    22 // we want to allocate the stats list,dlist arrays only once (or once per thread).
    23 // this function finds the largest array so we can allocate that max size when needed
     20// we want to allocate the StatDataSet arrays only once (or once per thread).  this
     21// function finds the largest value of Nmeasure so we can allocate that max size when
     22// needed
     23static int Nmax;
    2424void initMrel (Catalog *catalog, int Ncatalog) {
    2525
    … …  
    5454}
    5555
    56 void SetMrelInfoInit (SetMrelInfo *results, int allocLists) {
     56StatDataSet *StatDataSetAlloc (int Nsecfilt, int Nmax) {
     57
     58  int i;
     59
     60  StatDataSet *dataset = NULL;
     61  ALLOCATE (dataset, StatDataSet, Nsecfilt);
     62  for (i = 0; i < Nsecfilt; i++) {
     63    ALLOCATE (dataset[i].flxlist, double, Nmax);
     64    ALLOCATE (dataset[i].wgtlist, double, Nmax);
     65    ALLOCATE (dataset[i].errlist, double, Nmax);
     66    ALLOCATE (dataset[i].ranking,    int, Nmax);
     67    ALLOCATE (dataset[i].measSeq,    int, Nmax);
     68  } 
     69  return dataset;
     70}
     71
     72void StatDataSetFree (StatDataSet *dataset, int Nsecfilt) {
     73
     74  int i;
     75
     76  for (i = 0; i < Nsecfilt; i++) {
     77    FREE (dataset[i].flxlist);
     78    FREE (dataset[i].wgtlist);
     79    FREE (dataset[i].errlist);
     80    FREE (dataset[i].ranking);
     81    FREE (dataset[i].measSeq);
     82  } 
     83  FREE (dataset);
     84}
     85
     86void SetMrelInfoReset (SetMrelInfo *results) {
    5787  results->Nfew  = 0;
    5888  results->Ncode  = 0;
    … …  
    6292  results->Nmos  = 0;
    6393  results->Ngrid = 0;
     94}
     95
     96static int NsynthNames = 6;
     97static float synthMaxMags[] = {13.5, 13.5, 13.5, 13.0, 12.0, 13.5};
     98static char *synthNames[]   = { "g",  "r",  "i",  "z",  "y", "w"};
     99
     100void SetMrelInfoResetObject (SetMrelInfo *results) {
     101
     102  int i;
     103  for (i = 0; i < results->Nsecfilt; i++) {
     104    results->havePS1[i] = FALSE;
     105    results->haveSYN[i] = FALSE;
     106    results->needSYN[i] = FALSE;
     107    results->measSYN[i] = -1;
     108
     109    results->psfQfMax[i]  = NAN;
     110    results->psfQfPerfMax[i]  = NAN;
     111
     112    results->Nmeas[i] = 0;
     113    results->NmeasGood[i] = 0;
     114    results->Next[i] = 0;
     115    results->NexpPS1[i] = 0;
     116    results->havePS1[i] = 0;
     117    results->haveUbercal[i] = 0;
     118
     119    results->tessID[i]    = 0;
     120    results->projID[i]    = 0;
     121    results->skycellID[i] = 0;
     122
     123    results->minUbercalDist[i] = 1000;
     124
     125    results->psfData[i].Nlist = 0;
     126    results->aperData[i].Nlist = 0;
     127    results->kronData[i].Nlist = 0;
     128  }
     129}
     130
     131void SetMrelInfoInit (SetMrelInfo *results, int allocLists) {
     132
     133  int i;
     134  SetMrelInfoReset (results);
    64135  if (allocLists) {
    65     ALLOCATE (results->Mpsflist, double, Nmax);
    66     ALLOCATE (results->dpsflist, double, Nmax);
    67     ALLOCATE (results->wpsflist, double, Nmax);
     136    results->Nsecfilt = GetPhotcodeNsecfilt ();
     137
     138    results->psfData  = StatDataSetAlloc (results->Nsecfilt, Nmax);
     139    results->aperData = StatDataSetAlloc (results->Nsecfilt, Nmax);
     140    results->kronData = StatDataSetAlloc (results->Nsecfilt, Nmax);
     141
     142    ALLOCATE (results->psfqf_list,     double, Nmax);
     143    ALLOCATE (results->psfqfperf_list, double, Nmax);
     144    ALLOCATE (results->stargal_list,   double, Nmax);
     145
     146    ALLOCATE (results->havePS1,   int, results->Nsecfilt);
     147    ALLOCATE (results->haveSYN,   int, results->Nsecfilt);
     148    ALLOCATE (results->needSYN,   int, results->Nsecfilt);
     149    ALLOCATE (results->measSYN,   int, results->Nsecfilt);
     150    ALLOCATE (results->minSYN,  float, results->Nsecfilt);
     151
     152    ALLOCATE (results->psfQfMax,     float, results->Nsecfilt);
     153    ALLOCATE (results->psfQfPerfMax, float, results->Nsecfilt);
     154
     155    ALLOCATE (results->Nmeas,          int, results->Nsecfilt);
     156    ALLOCATE (results->NmeasGood,      int, results->Nsecfilt);
     157    ALLOCATE (results->Next,           int, results->Nsecfilt);
     158    ALLOCATE (results->NexpPS1,        int, results->Nsecfilt);
     159    ALLOCATE (results->haveUbercal,    int, results->Nsecfilt);
     160
     161    ALLOCATE (results->tessID,         int, results->Nsecfilt);
     162    ALLOCATE (results->projID,         int, results->Nsecfilt);
     163    ALLOCATE (results->skycellID,      int, results->Nsecfilt);
     164
     165    ALLOCATE (results->minUbercalDist, float, results->Nsecfilt);
     166
     167    SetMrelInfoResetObject (results);
     168
     169    for (i = 0; i < NsynthNames; i++) {
     170      int photcode = GetPhotcodeCodebyName (synthNames[i]);
     171      if (!photcode) continue;
     172      int Nsec = GetPhotcodeNsec (photcode);
     173      if (Nsec < 0) continue;
     174      results->minSYN[Nsec] = synthMaxMags[i];
     175    }
     176  } else {
     177    results->Nsecfilt = 0;
     178    results->psfData  = NULL;
     179    results->aperData = NULL;
     180    results->kronData = NULL;
     181
     182    results->psfqf_list     = NULL;
     183    results->psfqfperf_list = NULL;
     184    results->stargal_list   = NULL;
     185
     186    results->havePS1 = NULL;
     187    results->haveSYN = NULL;
     188    results->needSYN = NULL;
     189    results->measSYN = NULL;
     190    results->minSYN  = NULL;
     191
     192    results->psfQfMax = NULL;
     193    results->psfQfPerfMax = NULL;
     194
     195    results->Nmeas = NULL;
     196    results->NmeasGood = NULL;
     197    results->Next = NULL;
     198    results->NexpPS1 = NULL;
     199    results->haveUbercal = NULL;
     200
     201    results->tessID = NULL;
     202    results->projID = NULL;
     203    results->skycellID = NULL;
     204
     205    results->minUbercalDist = NULL;
    68206  }
    69207}
    70208
    71209void SetMrelInfoFree (SetMrelInfo *results) {
    72   free (results->Mpsflist);
    73   free (results->dpsflist);
    74   free (results->wpsflist);
     210  StatDataSetFree (results->psfData,  results->Nsecfilt);
     211  StatDataSetFree (results->aperData, results->Nsecfilt);
     212  StatDataSetFree (results->kronData, results->Nsecfilt);
     213  FREE (results->psfqf_list);
     214  FREE (results->psfqfperf_list);
     215  FREE (results->stargal_list);
     216
     217  FREE (results->havePS1);
     218  FREE (results->haveSYN);
     219  FREE (results->needSYN);
     220  FREE (results->measSYN);
     221  FREE (results->minSYN );
     222
     223  FREE (results->psfQfMax);
     224  FREE (results->psfQfPerfMax);
     225
     226  FREE (results->Nmeas);
     227  FREE (results->NmeasGood);
     228  FREE (results->Next);
     229  FREE (results->NexpPS1);
     230  FREE (results->haveUbercal);
     231
     232  FREE (results->tessID);
     233  FREE (results->projID);
     234  FREE (results->skycellID);
     235
     236  FREE (results->minUbercalDist);
    75237}
    76238
    … …  
    128290
    129291  for (i = 0; i < Ncatalog; i++) {
    130     setMrel_catalog_alt (catalog, i, -1, flatcorr, &results, Nsecfilt);
     292    setMrel_catalog_alt (catalog, i, FALSE, flatcorr, &results, Nsecfilt);
    131293    SetMrelInfoAccum (&summary, &results);
    132294  }
    … …  
    138300}
    139301
    140 int setMrelOutput (Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr) {
     302// setMrelOutput is NOT threaded because we read/write catalogs one at a time (for now)
     303int setMrelOutput (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
    141304
    142305  int i;
    143306
    144   // fprintf (stderr, "this version of setMrel is invalid for now\n");
    145   // exit (1);
    146        
    147307  int Nsecfilt = GetPhotcodeNsecfilt ();
    148308
    … …  
    151311  SetMrelInfoInit (&results, TRUE); // allocates results->list,dlist,wlist
    152312
    153   ALLOCATE (results.Maplist,  double, Nmax);
    154   ALLOCATE (results.daplist,  double, Nmax);
    155   ALLOCATE (results.waplist,  double, Nmax);
    156 
    157   ALLOCATE (results.Mkronlist, double, Nmax);
    158   ALLOCATE (results.dkronlist, double, Nmax);
    159   ALLOCATE (results.wkronlist, double, Nmax);
    160 
    161   ALLOCATE (results.psfqf_list, double, Nmax);
    162   ALLOCATE (results.psfqfperf_list, double, Nmax);
    163   ALLOCATE (results.stargal_list, double, Nmax);
    164 
    165313  for (i = 0; i < Ncatalog; i++) {
    166     setMrel_catalog_alt (catalog, i, pass, flatcorr, &results, Nsecfilt);
     314    setMrel_catalog_alt (catalog, i, TRUE, flatcorr, &results, Nsecfilt);
    167315    SetMrelInfoAccum (&summary, &results);
    168316  }
    … …  
    170318
    171319  SetMrelInfoFree (&results);
    172   free (results.Maplist);
    173   free (results.daplist);
    174   free (results.waplist);
    175 
    176   free (results.Mkronlist);
    177   free (results.dkronlist);
    178   free (results.wkronlist);
    179 
    180   free (results.psfqf_list);
    181   free (results.psfqfperf_list);
    182   free (results.stargal_list);
    183320  return (TRUE);
    184321}
    … …  
    261398
    262399  while (1) {
    263 
    264400    off_t i = getNextCatalogForThread(threadinfo->Ncatalog);
    265401    if (i == -1) {
    … …  
    271407    FlatCorrectionTable *flatcorr = threadinfo->flatcorr;
    272408
    273     // pass == -1 for anything other than the final pass
    274     setMrel_catalog_alt (catalog, i, -1, flatcorr, &results, Nsecfilt);
     409    setMrel_catalog_alt (catalog, i, FALSE, flatcorr, &results, Nsecfilt);
    275410    SetMrelInfoAccum (&threadinfo->summary, &results);
    276411  }
    … …  
    11421277            }
    11431278            Ndel ++;
     1279          } else {
     1280            catalog[i].measureT[ilist[k]].dbFlags &= ~ID_MEAS_POOR_PHOTOM;
     1281            if (final) {
     1282              // for the final pass, we have a duplicate set of values in measure and measureT
     1283              catalog[i].measure[ilist[k]].dbFlags &= ~ID_MEAS_POOR_PHOTOM;
     1284            }
    11441285          }
    11451286        }

  • trunk/Ohana/src/relphot/src/args.c

    r38441 r38466  
    324324  }
    325325
     326  // if -synthphot is chosen, an artificial w-band measurement is generated for each object
     327  // with r and i photometry.  this was used to apply the ubercal photometry to w-band images.
    326328  SyntheticPhotometry = FALSE;
    327329  if ((N = get_argument (argc, argv, "-synthphot"))) {
    … …  
    365367    mode = UPDATE_AVERAGES;
    366368  }
     369
     370  // if -synthphot-zpts is supplied, the map of corrections needed for the synthetic photometry is supplied
     371  // and used to tie down the synthetic magnitudes
     372  // NOTE: SYNTH_ZERO_POINTS is used for the same file in both -synthphot-zpts and -synthphot_means
    367373  SYNTH_ZERO_POINTS = NULL;
     374  if ((N = get_argument (argc, argv, "-synthphot-zpts"))) {
     375    remove_argument (N, &argc, argv);
     376    myAssert (N < argc, "missing argument to -synthphot-zpts");
     377    SYNTH_ZERO_POINTS = strcreate (argv[N]);
     378    remove_argument (N, &argc, argv);
     379  }
    368380  if ((N = get_argument (argc, argv, "-synthphot_means"))) {
    369381    mode = SYNTH_PHOT;
    … …  
    476488  }
    477489  SYNTH_ZERO_POINTS = NULL;
     490  if ((N = get_argument (argc, argv, "-synthphot-zpts"))) {
     491    remove_argument (N, &argc, argv);
     492    myAssert (N < argc, "missing argument to -synthphot-zpts");
     493    SYNTH_ZERO_POINTS = strcreate (argv[N]);
     494    remove_argument (N, &argc, argv);
     495  }
    478496  if ((N = get_argument (argc, argv, "-synthphot_means"))) {
    479497    MODE = MODE_SYNTH_PHOT;

  • trunk/Ohana/src/relphot/src/bcatalog.c

    r38441 r38466  
    3737  // flags used by the photometry analysis (excluding UBERCAL)
    3838  unsigned int PHOTOM_FLAGS =
    39     ID_MEAS_NOCAL |       // detection ignored for this analysis (photcode, time range)
    40     ID_MEAS_POOR_PHOTOM | // detection is photometry outlier
    41     ID_MEAS_SKIP_PHOTOM | // detection was ignored for photometry measurement
    42     ID_MEAS_AREA;         // detetion was outside acceptable area of device
    43 
     39    ID_MEAS_NOCAL          | // detection ignored for this analysis (photcode, time range)
     40    ID_MEAS_POOR_PHOTOM    | // detection is photometry outlier
     41    ID_MEAS_SKIP_PHOTOM    | // detection was ignored for photometry measurement
     42    ID_MEAS_AREA           | // detetion was outside acceptable area of device
     43    ID_MEAS_SYNTH_MAG      | // magnitude is synthetic
     44    ID_MEAS_STACK_PRIMARY  | // this stack measurement is in the primary skycell
     45    ID_MEAS_STACK_PHOT_SRC | // this measurement supplied the stack photometry
     46    ID_MEAS_PHOTOM_PSF     | // this measurement is used for the mean psf mag
     47    ID_MEAS_PHOTOM_APER    | // this measurement is used for the mean ap mag
     48    ID_MEAS_PHOTOM_KRON    ; // this measurement is used for the mean kron mag
     49   
    4450  /* exclude stars not in range or with too few measurements */
    4551  for (i = 0; i < catalog[0].Naverage; i++) {
    … …  
    5965      int Ns;
    6066
    61 # if (0)     
    62       DVOSecfiltFlags secfiltBits =
    63         ID_SECF_STAR_FEW |
    64         ID_SECF_STAR_POOR |
    65         ID_PHOTOM_PASS_0 |
    66         ID_PHOTOM_PASS_1 |
    67         ID_PHOTOM_PASS_2 |
    68         ID_PHOTOM_PASS_3 |
    69         ID_PHOTOM_PASS_4 |
    70         ID_SECF_USE_SYNTH |
    71         ID_SECF_USE_UBERCAL |
    72         ID_SECF_OBJ_EXT;
    73 # endif
    74 
    7567      for (Ns = 0; Ns < Nphotcodes; Ns++) {
    7668
    … …  
    7870        int Nsec = GetPhotcodeNsec(thisCode);
    7971
    80         // fprintf (stderr, "really use dvo_secfilt_init?");
    81         dvo_secfilt_init (&subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec]);
    82 
    83         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M      = NAN;
    84         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Map    = NAN;
    85         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].dM     = NAN;
    86         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Mstdev = NAN_S_SHORT;
    87         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Mchisq = NAN;
    88         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_20   = NAN_S_SHORT;
    89         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_80   = NAN_S_SHORT;
    90         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Ncode  = 0;
    91         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Nused  = 0;
    92         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~secfiltBits;
    93         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].ubercalDist = 1000;
     72        dvo_secfilt_init (&subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec], SECFILT_RESET_ALL);
    9473      }
    9574    }

  • trunk/Ohana/src/relphot/src/extra.c

    r38062 r38466  
    7070  return FALSE;
    7171}
     72
     73int is2MASS (int photcode) {
     74
     75  if ((photcode >= 2011) && (photcode <= 2013)) return TRUE;
     76  return FALSE;
     77}
     78
     79int isTYCHO (int photcode) {
     80
     81  if ((photcode == 2020) || (photcode == 2021)) return TRUE;
     82  return FALSE;
     83}
     84

  • trunk/Ohana/src/relphot/src/global_stats.c

    r35416 r38466  
    44
    55  StatType stN, stX, stS, imN, imX, imM, imD, msM, msX, msN, msD;
     6
     7  liststats_init (&stN);
     8  liststats_init (&stX);
     9  liststats_init (&stS);
     10  liststats_init (&imN);
     11  liststats_init (&imX);
     12  liststats_init (&imM);
     13  liststats_init (&imD);
     14  liststats_init (&msN);
     15  liststats_init (&msX);
     16  liststats_init (&msM);
     17  liststats_init (&msD);
    618
    719  // INITTIME;
    … …  
    2638 
    2739    fprintf (stderr, "   --- stats for %s ---\n", photcodes[Ns][0].name);
    28     fprintf (stderr, "meas / star:    %7.0f  %7.1f  %7.1f  %7.0f  %7.0f  %6d\n",   stN.median, stN.mean, stN.sigma, stN.min, stN.max, stN.Nmeas);
     40    fprintf (stderr, "meas / star:    %7.0f  %7.1f  %7.1f  %7.0f  %7.0f  %6d\n", stN.median, stN.mean, stN.sigma, stN.min, stN.max, stN.Nmeas);
    2941    fprintf (stderr, "dMrel star:     %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %6d\n", stS.median, stS.mean, stS.sigma, stS.min, stS.max, stS.Nmeas);
    30     fprintf (stderr, "chisq star:     %7.1f  %7.1f  %7.1f  %7.1f  %7.1f  %6d\n",   stX.median, stX.mean, stX.sigma, stX.min, stX.max, stX.Nmeas);
     42    fprintf (stderr, "chisq star:     %7.1f  %7.1f  %7.1f  %7.1f  %7.1f  %6d\n", stX.median, stX.mean, stX.sigma, stX.min, stX.max, stX.Nmeas);
    3143  }
    3244 

  • trunk/Ohana/src/relphot/src/liststats.c

    r34429 r38466  
    1616}
    1717
     18int liststats_init (StatType *stats) {
     19  stats->median  = NAN;
     20  stats->mean    = NAN;
     21  stats->sigma   = NAN;
     22  stats->error   = NAN;
     23  stats->chisq   = NAN;
     24  stats->min     = NAN;
     25  stats->max     = NAN;
     26  stats->Upper80 = NAN;
     27  stats->Lower20 = NAN;
     28  stats->total   = NAN;
     29  stats->Nmeas = 0;
     30  return TRUE;
     31}
     32
    1833int liststats (double *value, double *dvalue, double *weight, int N, StatType *stats) {
    1934 
    … …  
    2237
    2338  myAssert (stats->statmode != STATS_NONE, "programming error, liststats mode not set");
     39
     40  liststats_init (stats);
    2441
    2542  ke = ks = dMo = 0;

  • trunk/Ohana/src/relphot/src/plotstuff.c

    r35759 r38466  
    109109}
    110110
     111void plot_wait () {
     112  if (PLOTDELAY > 500000) {
     113    fprintf (stdout, "press return\n");
     114    if (fscanf (stdin, "%*c") != 1) fprintf (stderr, "\n");
     115  } else {
     116    usleep (PLOTDELAY);
     117  }
     118}
     119
    111120// plot the vector pair to a file with name defined by the varargs format
    112121void plot_list (Graphdata *graphdata, double *xlist, double *ylist, int N, char *label, char *format, ...) {
    … …  
    151160
    152161  if (SAVEPLOT) JpegPlot (graphdata, 0, filename);
    153   if (PLOTDELAY > 500000) {
    154     fprintf (stdout, "press return\n");
    155     if (fscanf (stdin, "%*c") != 1) fprintf (stderr, "\n");
    156   } else {
    157     usleep (PLOTDELAY);
    158   }
     162  plot_wait();
    159163  free (filename);
    160164}
    … …  
    166170  PlotVector (Npts, xlist, 0, 0, "x");
    167171  PlotVector (Npts, ylist, 1, 0, "y");
     172  plot_wait();
    168173}
    169174

  • trunk/Ohana/src/relphot/src/reload_catalogs.c

    r38441 r38466  
    3333  }
    3434
     35  // load the ZP corrections here
     36  if (SYNTH_ZERO_POINTS) SynthZeroPointsLoad (SYNTH_ZERO_POINTS);
     37
    3538  if (VERBOSE) fprintf (stderr, "re-loading catalog data\n");
    3639
    … …  
    6366
    6467    if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) {
    65       fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename);
    66       exit (1);
     68      fprintf (stderr, "ERROR: failure reading catalog %s, skipping\n", catalog.filename);
     69      continue;
     70      // exit (1);
    6771    }
    6872    if (VERBOSE && (catalog.Naverage_disk == 0)) {
    … …  
    168172    free (BOUNDARY_TREE);
    169173    BOUNDARY_TREE = tmppath;
     174  }
     175  if (SYNTH_ZERO_POINTS) {
     176    char *tmppath = abspath(SYNTH_ZERO_POINTS, DVO_MAX_PATH);
     177    free (SYNTH_ZERO_POINTS);
     178    SYNTH_ZERO_POINTS = tmppath;
    170179  }
    171180
    … …  
    234243    if (UPDATE_CATFORMAT) { strextend (&command, "-update-catformat %s", UPDATE_CATFORMAT); }
    235244    if (BOUNDARY_TREE)    { strextend (&command, "-boundary-tree %s", BOUNDARY_TREE); }
     245    if (SYNTH_ZERO_POINTS) { strextend {&command, "-synthphot-zpts %s", SYNTH_ZERO_POINTS); }
    236246
    237247    // deprecate

  • trunk/Ohana/src/relphot/src/relphot_objects.c

    r38441 r38466  
    3030  }
    3131
     32  // load the ZP corrections here
     33  if (SYNTH_ZERO_POINTS) SynthZeroPointsLoad (SYNTH_ZERO_POINTS);
     34
    3235  // load data from each region file, only use bright stars
    3336  for (i = 0; i < skylist[0].Nregions; i++) {
    … …  
    6063    if (RESET) {
    6164      Nsecfilt = catalog.Nsecfilt;
    62       DVOSecfiltFlags secfiltBits =
    63         ID_SECF_STAR_FEW    |
    64         ID_SECF_STAR_POOR   |
    65         ID_PHOTOM_PASS_0    |
    66         ID_PHOTOM_PASS_1    |
    67         ID_PHOTOM_PASS_2    |
    68         ID_PHOTOM_PASS_3    |
    69         ID_PHOTOM_PASS_4    |
    70         ID_SECF_USE_SYNTH   |
    71         ID_SECF_USE_UBERCAL |
    72         ID_SECF_HAS_STACK   |
    73         ID_SECF_HAS_PS1     |
    74         ID_SECF_OBJ_EXT     ;
    75      
     65
    7666      DVOAverageFlags astromBits =
    7767        ID_STAR_FIT_AVE   |
    … …  
    8979        catalog.average[j].stargal   = NAN;     // reset (will be re-calculated here)
    9080        for (k = 0; k < Nsecfilt; k++) {
    91           dvo_secfilt_init (&catalog.secfilt[j*Nsecfilt + k]);
    92           catalog.secfilt[j*Nsecfilt + k].flags &= ~secfiltBits; // XXX: only reset photom bits:
     81          dvo_secfilt_init (&catalog.secfilt[j*Nsecfilt + k], SECFILT_RESET_ALL);
    9382        }
    9483      }
    … …  
    9988    initMrel (&catalog, 1);
    10089    setMrelFinal (&catalog, NULL, TRUE);
     90    // XXX if we want to have options for setting warp, chip, stack independently, we need to init only the desired ones
    10191
    10292    if (!UPDATE) {
    … …  
    164154    if (UPDATE)        { strextend (&command, "-update"); }
    165155    if (!KEEP_UBERCAL) { strextend (&command, "-reset-ubercal"); }
     156    if (SYNTH_ZERO_POINTS) { snprintf (tmpline, 1024, "%s -synthphot-zpts %s", command, SYNTH_ZERO_POINTS); strcpy (command, tmpline); }
    166157
    167158    // if (SET_MREL_VERSION != 1) { strextend (&command, "-set-mrel-version %d", SET_MREL_VERSION); }

  • trunk/Ohana/src/relphot/src/relphot_synthphot.c

    r38441 r38466  
    2828
    2929  // load the ZP corrections here
    30   SynthZeroPoints *zpts = SynthZeroPointsLoad (SYNTH_ZERO_POINTS);
     30  SynthZeroPointsLoad (SYNTH_ZERO_POINTS);
     31  SynthZeroPoints *zpts = SynthZeroPointsGet ();
    3132  if (!zpts) exit (2);
    3233

  • trunk/Ohana/src/relphot/src/setMrelCatalog.c

    r38153 r38466  
    11# include "relphot.h"
     2int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags flags);
     3
    24# define UBERCAL_WEIGHT 100.0
    3 
    4 # define SKIP_THIS_MEAS(REASON) {                               \
    5     measureT[k].dbFlags |= ID_MEAS_SKIP_PHOTOM; \
    6     if (measure) {                                      \
    7       measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM;                \
    8     }                                                           \
    9     results->REASON ++;                                         \
     5void sort_StatDataSet (StatDataSet *dataset);
     6
     7# define SKIP_THIS_MEAS(REASON) {                       \
     8    results->REASON ++;                                 \
    109    continue; }
    1110
    1211# define SKIP_THIS_MEAS_STACK(REASON) {                 \
    13     measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM;          \
    1412    continue; }
     13
     14# define CHECK_VALID_MAG(MAG,D_MAG) (isfinite(MAG) && isfinite(D_MAG) && (MAG > -5.0) && (MAG < 30.0))
     15# define CHECK_VALID_FLUX(FLUX,D_FLUX) (isfinite(FLUX) && isfinite(D_FLUX))
    1516
    1617static float MagToFlux (float Mag) {
    … …  
    3637}
    3738
    38 int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
     39int setMrel_catalog_alt (Catalog *catalog, int Nc, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
    3940
    4041  off_t j;
    4142
    42   liststats_setmode (&results->psfstats, STATMODE);
    43   liststats_setmode (&results->apstats, STATMODE);
     43  liststats_setmode (&results->psfstats,  STATMODE);
     44  liststats_setmode (&results->apstats,   STATMODE);
    4445  liststats_setmode (&results->kronstats, STATMODE);
    4546
    46   SetMrelInfoInit (results, FALSE); // do not allocate list,dlist,wlist arrays
    47 
    48   int isSetMrelFinal = (pass >= 0);
    49 
    50   // XX char *primaryCell = NULL;
    51   // XX if (isSetMrelFinal) {
    52   // XX   ALLOCATE (primaryCell, char, DVO_MAX_PATH);
    53   // XX }
     47  SetMrelInfoReset (results); // reset the countesrs
    5448
    5549  for (j = 0; j < catalog[Nc].Naverage; j++) {
    56     // we are guaranteed to have averageT, but not average
    57     off_t m = catalog[Nc].averageT[j].measureOffset;
    58     Average *average = catalog[Nc].average ? &catalog[Nc].average[j] : NULL;
    59     Measure *measure = catalog[Nc].measure ? &catalog[Nc].measure[m] : NULL;
    60     SecFilt *secfilt = catalog[Nc].secfilt ? &catalog[Nc].secfilt[j*Nsecfilt] : NULL;
    61     setMrelAverageExposure (m, Nc, pass, flatcorr, results, average, &catalog[Nc].averageT[j], secfilt, measure, &catalog[Nc].measureT[m], &catalog[Nc].found_t[Nsecfilt*j]);
     50
     51    setMrelAverageExposure (&catalog[Nc], Nc, j, Nsecfilt, isSetMrelFinal, flatcorr, results);
    6252
    6353    // only apply Stack operation on setMrelFinal in first pass
    64     if (isSetMrelFinal && (pass == 0) && !IS_DIFF_DB) {
    65       setMrelAverageStack (m, Nc, flatcorr, results, average, secfilt, &catalog[Nc].measure[m], catalog[Nc].Nmeasure);
    66       setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]);
     54    if (isSetMrelFinal && !IS_DIFF_DB) {
     55      setMrelAverageStack (&catalog[Nc], Nc, j, Nsecfilt, flatcorr);
    6756    }
    6857
    6958    // only measure force-warp mean values if issetMrelFinal (make it optional?)
    7059    if (isSetMrelFinal) {
    71       setMrelAverageForcedWarp (m, Nc, pass, results, average, secfilt, measure, &catalog[Nc].foundWarp_t[Nsecfilt*j], catalog[Nc].Nmeasure);
    72       setGlobalObjStats (average, measure);
    73     }
    74   }
    75   // if (primaryCell) free (primaryCell);
     60      setMrelAverageForcedWarp (&catalog[Nc], Nc, j, Nsecfilt, flatcorr, results);
     61    }
     62  }
     63
    7664  return (TRUE);
    7765}
    7866
     67// NOTE:
     68// Msys is measure[i].M + zp corrections
     69// Mcal is image[j].Mcal
     70// Mmos and Mgrid are offsets for mosaic and grid
     71
     72// tie down reference photometry if the -refcode (code) option is selected
     73// eg, -refcode g_SDSS
     74// this probably makes no sense in the context of multifilter analysis
     75// XXX probably need to use the photocde table to assign reference mag weights.
     76
     77// dlist gives the error per measurement, wlist gives the weight
     78// we can modify the error and weight in a few ways:
     79// 1) MIN_ERROR guarantees a floor
     80// 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
     81// 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
     82// 4) some reference photcode of some kind can be specified as fixed and have a high weight
     83
    7984// set mean of chip measurements (selected by photcode range for now):
    80 int setMrelAverageExposure (off_t measureOffset, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found) {
    81 
     85int setMrelAverageExposure (Catalog *catalog, int cat, off_t ave, int Nsecfilt, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results) {
     86
     87  int Nsec;
    8288  off_t k;
    83   float Msys = 0, Mcal= 0, Mmos = 0, Mgrid = 0;
    84 
    85   // we are measuring means for 3 types of mags: psf, ap, kron.  I am using the psf mag
    86   // error for the ap mags, but krons have their own errors.  it is an open question if I
    87   // should be doing weighted or unweighted fits (this is a user option)
    88   double *Mpsflist  = results->Mpsflist;
    89   double *dpsflist  = results->dpsflist;
    90   double *wpsflist  = results->wpsflist;
    91 
    92   double *Maplist   = results->Maplist;
    93   double *daplist   = results->daplist;
    94   double *waplist   = results->waplist;
    95 
    96   double *Mkronlist = results->Mkronlist;
    97   double *dkronlist = results->dkronlist;
    98   double *wkronlist = results->wkronlist;
    99 
    100   StatType *psfstats  = &results->psfstats;
    101   StatType *apstats   = &results->apstats;
    102   StatType *kronstats = &results->kronstats;
    103 
    104   int isSetMrelFinal = (pass >= 0);
     89
     90  // we are guaranteed to have averageT, but not average
     91  AverageTiny *averageT = &catalog[0].averageT[ave];
     92
     93  off_t measureOffset = averageT->measureOffset;
     94  MeasureTiny *measureT = &catalog[0].measureT[measureOffset];
     95
     96  // we are NOT guaranteed to have average, measure, secfilt
     97  Average *average  = catalog[0].average     ? &catalog[0].average[ave]               : NULL;
     98  Measure *measure  = catalog[0].measure     ? &catalog[0].measure[measureOffset]     : NULL;
     99  SecFilt *secfilt  = catalog[0].secfilt     ? &catalog[0].secfilt[ave*Nsecfilt]      : NULL;
     100  char *measureRank = catalog[0].measureRank ? &catalog[0].measureRank[measureOffset] : NULL;
     101
     102  // we are measuring means for 3 types of mags: psf, ap, kron.
     103
     104  // in the final assignment, set the mean mags even if only 1 measurements exists
     105  int Nminmeas = isSetMrelFinal ? 1 : STAR_TOOFEW + 1;
     106
     107  // isSetMrelFinal : in the final pass, we set psf, kron, ap mag values and extra stats; other passes only do psf mags
    105108
    106109  // option for a test print
    107   if (FALSE && (average[0].objID == 0x7146) && (average[0].catID == 0x49d8)) {
     110  if (FALSE && (averageT[0].objID == 0x7146) && (averageT[0].catID == 0x49d8)) {
    108111    fprintf (stderr, "test obj\n");
    109112    print_measure_set_alt (average, secfilt, measure);
    110113  }
    111114
     115  SetMrelInfoResetObject (results); // reset the per-object arrays
     116
     117  int NextPS1 = 0;
     118  int NpsfPS1 = 0;
    112119  int GoodPS1 = FALSE;
    113120  int Good2MASS = FALSE;
    114121  int Galaxy2MASS = FALSE;
    115 
    116   int NextPS1 = 0;
    117   int NpsfPS1 = 0;
    118 
    119   int Ns;
    120   for (Ns = 0; Ns < Nphotcodes; Ns++) {
    121 
    122     int thisCode = photcodes[Ns][0].code;
    123     int Nsec = GetPhotcodeNsec(thisCode);
    124 
    125     /* calculate the average mag in this SEC photcode for a single star */
    126 
    127     /* star/photcodes already calibrated */
    128     if ( isSetMrelFinal && found[Nsec]) continue; 
    129      
    130     // skip bad stars
    131     if (!isSetMrelFinal && (secfilt[Nsec].flags & STAR_BAD)) continue;
    132 
    133     int NexpPS1 = 0;
    134     int Ncode = 0;
    135     int Next = 0;
    136     int haveSynth = FALSE;
    137     int havePS1   = FALSE;
    138 
    139     int forceSynth = FALSE;
    140     int forceSynthEntry = -1;
    141 
    142     int haveUbercal = FALSE;
    143 
    144     int minUbercalDist = 1000;
     122  int haveTYCHO = FALSE;
     123
     124  float stargalmax = 0.0;
     125
     126  // assign measurements to the photcode lists
     127  for (k = 0; k < averageT[0].Nmeasure; k++) {
     128
     129    // these bits should not be set unless we use them in this pass
     130    measure[k].dbFlags &= ~(ID_MEAS_PHOTOM_PSF | ID_MEAS_PHOTOM_APER | ID_MEAS_PHOTOM_KRON);
     131
     132    // skip measurements that do not have a valid photcode (raise exception?)
     133    PhotCode *code = GetPhotcodebyCode (measureT[k].photcode);
     134    myAssert (code, "invalid photcode??");
     135
     136    // SKIP gpc1 stack data (hard-wired photcodes)
     137    if (isGPC1stack(measureT[k].photcode)) continue;
    145138   
    146     off_t meas = measureOffset;
    147 
    148     float psfQfMax = 0.0;
    149     float psfQfPerfMax = 0.0;
    150 
    151     int Nap = 0;
    152     int Npsf = 0;
    153     int Nkron = 0;
    154     for (k = 0; k < averageT[0].Nmeasure; k++, meas++) {
    155 
    156       // skip measurements that do not match the current photcode
    157       PhotCode *code = GetPhotcodebyCode (measureT[k].photcode);
    158       if (!code) continue;
    159       if (code->equiv != thisCode) { continue; }
    160       Ncode ++;
    161 
    162       // are we a PS1 exposure photcode?
    163       if (isGPC1chip(measureT[k].photcode)) NexpPS1 ++;
    164 
    165       // SKIP gpc1 stack data
    166       if (isGPC1stack(measureT[k].photcode)) continue;
    167 
    168       // SKIP gpc1 forced-warp data
    169       if (isGPC1warp(measureT[k].photcode)) continue;
    170 
    171       if (isSetMrelFinal) {
    172         if (measure[k].psfQF     > psfQfMax)     psfQfMax     = measure[k].psfQF;
    173         if (measure[k].psfQFperf > psfQfPerfMax) psfQfPerfMax = measure[k].psfQFperf;
    174       }
    175 
    176       if (measureT[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad);
    177 
    178       if (getImageEntry (meas, cat) < 0) {
    179         // measurements without an image are either external reference photometry or
    180         // data for which the associated image has not been loaded (probably because of
    181         // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
    182         Mmos = Mgrid = 0;
    183         Mcal = measureT[k].Mcal; // check that this is zero for loaded REF value
     139    // SKIP gpc1 forced-warp data (hard-wired photcodes)
     140    if (isGPC1warp(measureT[k].photcode)) continue;
     141   
     142    if (isTYCHO(measureT[k].photcode)) { haveTYCHO = TRUE; }
     143
     144    if (is2MASS(measureT[k].photcode)) {
     145      if (measureT[k].photFlags & 0x00c00000) {
     146        Galaxy2MASS = TRUE; // per object value
     147      }
     148      if (measureT[k].photFlags & 0x00000007) {
     149        Good2MASS = TRUE; // per object value
     150      }
     151    }
     152
     153    // assign the Nsec value so we can assign to the right lists
     154    int Nsec = GetPhotcodeNsec (code->equiv);
     155    if (Nsec < 0) continue; // skip measurements which do not have an equiv average photcode
     156
     157    results->Nmeas[Nsec] ++;
     158
     159    // various steps only relevant to the final pass
     160    if (isSetMrelFinal) {
     161      if (isfinite(measure[k].psfQF)     && (measure[k].psfQF     > results->psfQfMax[Nsec]))     results->psfQfMax[Nsec]     = measure[k].psfQF;
     162      if (isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > results->psfQfPerfMax[Nsec])) results->psfQfPerfMax[Nsec] = measure[k].psfQFperf;
     163
     164      // are we a PS1 exposure photcode? (hard-wired photcodes)
     165      if (isGPC1chip(measure[k].photcode)) {
     166        results->NexpPS1[Nsec] ++;
     167        results->havePS1[Nsec] = TRUE;
     168       
     169        // only count psfQF > 0.85 here
     170        if (isfinite(measure[k].extNsigma) && isfinite(measure[k].psfQF) && (measure[k].psfQF > 0.85)) {
     171          stargalmax = MAX (stargalmax, measure[k].extNsigma);
     172        }
     173      }
     174
     175      // force the use of SYN even if we have PS1 mags?
     176      if (isGPC1synth(measure[k].photcode)) {
     177        results->haveSYN[Nsec] = TRUE;
     178        results->measSYN[Nsec] = k;
     179        measureT[k].dbFlags |= ID_MEAS_SYNTH_MAG; // redundant with photcode
     180        // note that synthetic mags have the real mags in measure.M (no zero point offset)
     181        if (measure[k].M < results->minSYN[Nsec]) results->needSYN[Nsec] = TRUE;
     182        continue;
     183      }
     184    }
     185
     186    // ** Choose the calibration (depends on the mode : do I have an image reference or not?)
     187    float Mcal = 0, Mmos = 0, Mgrid = 0;
     188    off_t meas = measureOffset + k;
     189    if (getImageEntry (meas, cat) < 0) {
     190      // measurements without an image are either external reference photometry or
     191      // data for which the associated image has not been loaded (probably because of
     192      // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
     193      Mmos = Mgrid = 0;
     194      Mcal = measureT[k].Mcal; // check that this is zero for loaded REF value
     195    } else {
     196      // getMcal_alt returns image[].Mcal modified by flatcorr(image.photom_map_id,x,y)
     197      Mcal  = getMcal_alt (meas, cat, flatcorr, measureT[k].Xccd, measureT[k].Yccd);
     198      if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
     199      Mmos  = getMmos  (meas, cat);
     200      if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
     201      Mgrid = getMgrid (meas, cat);
     202      if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
     203    }
     204
     205    int myUbercalDist = getUbercalDist(meas, cat);
     206    results->minUbercalDist[Nsec] = MIN(results->minUbercalDist[Nsec], myUbercalDist);
     207
     208    int isUbercal = (measureT[k].dbFlags & ID_MEAS_PHOTOM_UBERCAL);
     209
     210    if (isUbercal) results->haveUbercal[Nsec] = TRUE; // haveUbercal is set per secfilt, isUbercal is per measure XXX define this array
     211
     212    int useUbercalWeight = isUbercal || (refPhotcode && (code->code == refPhotcode->code));
     213
     214    float Map  = NAN;
     215    float dMap = NAN;
     216    if (isSetMrelFinal) {
     217      Map = PhotCat (&measure[k], MAG_CLASS_APER);
     218      dMap = MAX (hypot(measure[k].dMap,   code->photomErrSys), MIN_ERROR);
     219      if (CHECK_VALID_MAG(Map, dMap)) {
     220        int Nap = results->aperData[Nsec].Nlist;
     221        results->aperData[Nsec].flxlist[Nap] = Map - Mcal - Mmos - Mgrid;
     222        results->aperData[Nsec].errlist[Nap] = dMap;
     223        results->aperData[Nsec].wgtlist[Nap] = useUbercalWeight ? UBERCAL_WEIGHT : 1.0;
     224        results->aperData[Nsec].ranking[Nap] = measureRank[k];
     225        results->aperData[Nsec].measSeq[Nap] = k;
     226        results->aperData[Nsec].Nlist ++;
     227      }
     228    }
     229
     230    float Mkron, dMkron;
     231    if (isSetMrelFinal) {
     232      Mkron = PhotCat (&measure[k], MAG_CLASS_KRON);
     233      dMkron = MAX (hypot(measure[k].dMkron, code->photomErrSys), MIN_ERROR);
     234      if (CHECK_VALID_MAG(Mkron, dMkron)) {
     235        int Nkron = results->kronData[Nsec].Nlist;
     236        results->kronData[Nsec].flxlist[Nkron] = Mkron - Mcal - Mmos - Mgrid;
     237        results->kronData[Nsec].errlist[Nkron] = dMkron;
     238        results->kronData[Nsec].wgtlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0;
     239        results->kronData[Nsec].ranking[Nkron] = measureRank[k];
     240        results->kronData[Nsec].measSeq[Nkron] = k;
     241        results->kronData[Nsec].Nlist ++;
     242      }
     243    }
     244
     245    float Mpsf = PhotSysTiny (&measureT[k], &averageT[0], &secfilt[0], MAG_CLASS_PSF);
     246    float dMpsf  = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR);
     247    if (CHECK_VALID_MAG(Mpsf, dMpsf)) {
     248      int Npsf = results->psfData[Nsec].Nlist;
     249      results->psfData[Nsec].flxlist[Npsf] = Mpsf - Mcal - Mmos - Mgrid;
     250      results->psfData[Nsec].errlist[Npsf] = dMpsf;
     251      results->psfData[Nsec].wgtlist[Npsf] = isUbercal ? UBERCAL_WEIGHT : 1.0;
     252      results->psfData[Nsec].ranking[Npsf] = measureRank ? measureRank[k] : 0;
     253      results->psfData[Nsec].measSeq[Npsf] = k;
     254      results->psfData[Nsec].Nlist ++;
     255    }
     256
     257    // count the extended detections (all PS1 bands)
     258    if (isSetMrelFinal && isGPC1chip(measureT[k].photcode) && !isnan(Map) && !isnan(Mpsf)) {
     259      float dMagAp = Mpsf - Map;
     260      float SigmaAp = hypot(0.1, 2.5*dMpsf);
     261      // XXX this is still quite ad hoc, but at least it:
     262      // (a) converges to 0.1 mag offset at the bright end
     263      // (b) converges to 0.5 mag offset at the faint end (dM = 0.2)
     264      if (dMagAp > SigmaAp) {
     265        results->Next[Nsec] ++;
     266        NextPS1 ++;
    184267      } else {
    185         Mcal  = getMcal_alt (meas, cat, flatcorr, measureT[k].Xccd, measureT[k].Yccd);
    186         if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
    187         Mmos  = getMmos  (meas, cat);
    188         if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
    189         Mgrid = getMgrid (meas, cat);
    190         if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
    191       }
    192 
    193       // skip some absurd values NAN, < 0.0, > 30.0
    194       Msys = PhotSysTiny (&measureT[k], &averageT[0], &secfilt[0], MAG_CLASS_PSF);
    195       if (isnan(Msys)) SKIP_THIS_MEAS(Nsys);
    196       if (Msys <  0.0) SKIP_THIS_MEAS(Nsys);
    197       if (Msys > 30.0) SKIP_THIS_MEAS(Nsys);
    198 
    199       int myUbercalDist = getUbercalDist(meas, cat);
    200       minUbercalDist = MIN(minUbercalDist, myUbercalDist);
    201 
    202       int isUbercal = (measureT[k].dbFlags & ID_MEAS_PHOTOM_UBERCAL);
    203       float dMpsf = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR);
    204 
    205       if (isUbercal) haveUbercal = TRUE; // haveUbercal is set per secfilt, isUbercal is per measure
    206 
    207       if (isSetMrelFinal) {
    208         // special options for PS1 data
    209         if ((measure[k].photcode >= 10000) && (measure[k].photcode <= 10500)) {
    210           // count the extended detections
    211           if (!isnan(measure[k].Map)) {
    212             float dMagAp = measure[k].M - measure[k].Map;
    213             float SigmaAp = hypot(0.1, 2.5*measure[k].dM);
    214             // XXX this is still quite ad hoc, but at least it:
    215             // (a) converges to 0.1 mag offset at the bright end
    216             // (b) converges to 0.5 mag offset at the faint end (dM = 0.2)
    217             if (dMagAp > SigmaAp) {
    218               Next ++;
    219               NextPS1 ++;
    220             } else {
    221               NpsfPS1 ++;
    222             }
    223           }
    224           havePS1 = TRUE;
     268        NpsfPS1 ++;
     269      }
     270    }
     271  }
     272
     273  SynthZeroPoints *synthzpts = SynthZeroPointsGet();
     274
     275  // find max values across filters
     276  float psfQfMax     = 0.0;
     277  float psfQfPerfMax = 0.0;
     278
     279  // now calculate the mean stats for the Nsec bands.
     280  for (Nsec = 0; Nsec < Nsecfilt; Nsec++) {
     281    dvo_secfilt_init (&secfilt[Nsec], SECFILT_RESET_CHIP); // this does not reset astrometry or STACK bits
     282
     283    // XXX hardwired
     284    if ((Nsec < 4) || (Nsec == 8)) {
     285      secfilt[Nsec].Ncode = results->NexpPS1[Nsec];
     286    } else {
     287      secfilt[Nsec].Ncode = results->Nmeas[Nsec]; // 2MASS data if it exists
     288    }
     289
     290    if (haveTYCHO) {
     291      secfilt[Nsec].flags |= ID_SECF_HAS_TYCHO;
     292    }
     293
     294    // force the synthetic magnitude (apply zpt calibration if no tycho photometry exists)
     295    if (results->needSYN[Nsec]) {
     296      if (!isSetMrelFinal) continue;
     297
     298      // use the single SYNTH value instead of the other mags here
     299      float Mpsf = measure[results->measSYN[Nsec]].M;
     300      float ZP = 0.0;
     301
     302      // if we have loaded the synth zpt correction, and do NOT have a tycho measurement, apply the correction
     303      if (!haveTYCHO && synthzpts) {
     304        // need to look up the (X,Y) coords from this (R,D) location
     305        double X, Y;
     306        double R = ohana_normalize_angle_to_midpoint (measure->R, 0.0); // XXX 0.0 or 180.0 for center?
     307        RD_to_XY (&X, &Y, R, measure->D, &synthzpts->coords);
     308        if (X < 0) return FALSE;
     309        if (Y < 0) return FALSE;
     310        if (X >= synthzpts->Nx) return FALSE;
     311        if (Y >= synthzpts->Ny) return FALSE;
     312       
     313        int Xpix = X;
     314        int Ypix = Y;
     315        int Npix = Xpix + Ypix*synthzpts->Nx;
     316       
     317        float *value = (float *) synthzpts->matrix[Nsec].buffer;
     318        ZP = !isnan(value[Npix]) ? value[Npix] : 0.0;
     319        if (isfinite(value[Npix])) {
     320          secfilt[Nsec].flags |= ID_SECF_FIX_SYNTH;
    225321        }
    226 
    227         // count extended detections for 2MASS (XXX NOTE hardwired photcodes 2011, 2012, 2013)
    228         if ((measure[k].photcode >= 2011) && (measure[k].photcode <= 2013)) {
    229           if (measure[k].photFlags & 0x00c00000) {
    230             Next ++;
    231             Galaxy2MASS = TRUE;
    232           }
    233           if (pass == 0) {
    234             if (measure[k].photFlags & 0x00000007) {
    235               Good2MASS = TRUE;
    236             } else {
    237               // detections without one of these bits should only be used in PASS_1
    238               SKIP_THIS_MEAS(Nbad);
    239             }
    240           }
    241         }
    242 
    243         // Blindly accepth the SYNTH mags if we are above saturation, otherwise,
    244         // ignore SYNTH photcodes until PASS == 4 (where we also accept saturated stars)
    245         if ((measure[k].photcode >= 3001) && (measure[k].photcode <= 3005)) {
    246           // something of a hack: force object to use synth values if synth mags >>
    247           // saturation (3pi instrumental mags < -15)
    248           float MaxMagForceSynth = NAN;
    249           switch (measure[k].photcode) {
    250             case 3001:
    251               MaxMagForceSynth = 13.64;
    252               break;
    253             case 3002:
    254               MaxMagForceSynth = 13.76;
    255               break;
    256             case 3003:
    257               MaxMagForceSynth = 13.74;
    258               break;
    259             case 3004:
    260               MaxMagForceSynth = 12.94;
    261               break;
    262             case 3005:
    263               MaxMagForceSynth = 12.01;
    264               break;
    265           }
    266           if (measureT[k].M < MaxMagForceSynth) {
    267             forceSynth = TRUE;
    268             forceSynthEntry = Npsf;
    269           } else {
    270             if (pass < 4) {
    271               SKIP_THIS_MEAS(Nbad);
    272             }
    273             haveSynth = TRUE;
    274           }
    275         }
    276 
    277         // Map (Maplist) and Mkron (Mkronlist,dkronlist) are used to calculate mean mags per filter
    278         float Map = PhotCat (&measure[k], MAG_CLASS_APER);
    279         if (!isnan(Map)) {
    280             Maplist[Nap] = Map - Mcal - Mmos - Mgrid;
    281             daplist[Nap] = dMpsf;// XXX check on this...
    282             waplist[Nap] = isUbercal ? UBERCAL_WEIGHT : 1.0;
    283             Nap ++;
    284         }
    285 
    286         float Mkron = PhotCat (&measure[k], MAG_CLASS_KRON);
    287         if (!isnan(Mkron)) {
    288           Mkronlist[Nkron] = Mkron - Mcal - Mmos - Mgrid;
    289           dkronlist[Nkron] = measure[k].dMkron;
    290           wkronlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0;
    291           Nkron ++;
    292         }
    293       } // if (isSetMrelFinal)
    294 
    295       // dlist gives the error per measurement, wlist gives the weight
    296       // we can modify the error and weight in a few ways:
    297       // 1) MIN_ERROR guarantees a floor
    298       // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
    299       // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
    300       // 4) some reference photcode of some kind can be specified as fixed and have a high weight
    301       Mpsflist[Npsf] = Msys - Mcal - Mmos - Mgrid;
    302       dpsflist[Npsf] = dMpsf;
    303       wpsflist[Npsf] = isUbercal ? UBERCAL_WEIGHT : 1.0;
    304 
    305       // NOTE:
    306       // Msys is measure[i].M + zp corrections
    307       // Mcal is image[j].Mcal
    308       // Mmos and Mgrid are offsets for mosaic and grid
    309 
    310       // tie down reference photometry if the -refcode (code) option is selected
    311       // eg, -refcode g_SDSS
    312       // this probably makes no sense in the context of multifilter analysis
    313       // XXX probably need to use the photocde table to assign reference mag weights.
    314       if (refPhotcode) {
    315         if (code->code == refPhotcode->code) {
    316           wpsflist[Npsf] = UBERCAL_WEIGHT;
    317         }
    318       }
    319       Npsf ++;
    320     }
    321 
    322     int Nminmeas = isSetMrelFinal ? 1 : STAR_TOOFEW + 1;
    323 
    324     // XXX : ugh : another hard-wired photcode entry...
    325     if (isSetMrelFinal && (pass == 0)) {
    326       if ((thisCode < 6) || (thisCode == 9)) {
    327         secfilt[Nsec].Ncode = NexpPS1;
    328       } else {
    329         secfilt[Nsec].Ncode = Ncode; // 2MASS data if it exists
    330       }
    331     }
    332 
    333     // when performing the grid analysis, STAR_TOOFEW should be set to 1;
    334     if (Npsf < Nminmeas) { /* too few measurements */
    335       // fprintf (f, "%10.6f %10.6f %d %d %d\n", averageT[0].R, averageT[0].D, measureT[0].imageID, Npsf, STAR_TOOFEW);
     322      }
     323      secfilt[Nsec].M      = Mpsf + ZP;
     324      secfilt[Nsec].dM     = 0.6;
     325      secfilt[Nsec].Mchisq = 0.0;
     326      secfilt[Nsec].flags |= ID_SECF_USE_SYNTH;
     327      continue;
     328    }
     329
     330    // if too few valid measurements meet the minimum criteria, go to the next entry
     331    StatType *psfstats = &results->psfstats;
     332    int Nranking = magStatsByRanking (&results->psfData[Nsec], psfstats);
     333    if (Nranking < Nminmeas) {
    336334      secfilt[Nsec].flags |= ID_STAR_FEW;
    337       if (Ncode == 0) {
    338         results->Ncode ++;
    339       } else {
    340         results->Nfew ++;
    341       }
    342       continue;
    343335    } else {
    344       secfilt[Nsec].flags &= ~ID_STAR_FEW;
    345     }   
    346 
    347     if (forceSynth) {
    348       // use the single SYNTH value instead of the other mags here
    349       myAssert ((forceSynthEntry < Npsf) && (forceSynthEntry >= 0), "programming error");
    350       Mpsflist[0] = Mpsflist[forceSynthEntry];
    351       dpsflist[0] = dpsflist[forceSynthEntry];
    352       wpsflist[0] = wpsflist[forceSynthEntry];
    353       Npsf = 1;
    354     }
    355     liststats (Mpsflist, dpsflist, wpsflist, Npsf, psfstats);
    356 
    357     secfilt[Nsec].M      = psfstats->mean;
    358     secfilt[Nsec].dM     = psfstats->error;
    359     secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
    360 
    361     // when running -averages, we have no information about the images, so we cannot set this
    362     if (minUbercalDist > -1) {
    363       secfilt[Nsec].ubercalDist = minUbercalDist;
    364     }
     336      secfilt[Nsec].M    = psfstats->mean;
     337      secfilt[Nsec].dM   = psfstats->error;
     338      secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
     339    }
     340    int minRank = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10;
    365341
    366342    if (isSetMrelFinal) {
    367       found[Nsec] = TRUE;
    368 
    369       secfilt[Nsec].Mstdev = psfstats->sigma; // Mstdev is in millimags (not enough space for more precision)
    370       // secfilt[Nsec].Ncode = Ncode;
    371       secfilt[Nsec].Nused = psfstats->Nmeas;
    372 
    373       secfilt[Nsec].Mmax = psfstats->max;
    374       secfilt[Nsec].Mmin = psfstats->min;
    375 
    376       secfilt[Nsec].psfQfMax     = psfQfMax;
    377       secfilt[Nsec].psfQfPerfMax = psfQfPerfMax;
    378 
    379       // NOTE : use the modified weight for apmags as well as psf mags
    380       liststats (Maplist, daplist, waplist, Nap, apstats);
    381       secfilt[Nsec].Map  = Nap > 0 ? apstats->mean : NAN;
    382       secfilt[Nsec].dMap  = Nap > 0 ? apstats->error : NAN;
    383       secfilt[Nsec].sMap  = Nap > 0 ? apstats->sigma : NAN;
    384       secfilt[Nsec].NusedAp  = Nap;
    385 
    386       liststats (Mkronlist, dkronlist, wkronlist, Nkron, kronstats);
    387       secfilt[Nsec].Mkron  = Nkron > 0 ? kronstats->mean  : NAN;
    388       secfilt[Nsec].dMkron = Nkron > 0 ? kronstats->error : NAN;
    389       secfilt[Nsec].sMkron = Nkron > 0 ? kronstats->sigma : NAN;
    390       secfilt[Nsec].NusedKron  = Nkron;
    391 
    392       // NOTE: for 2MASS measurements, Next should be 1, as should N
    393       if ((Next > 0) && (Next > 0.5*Npsf)) {
     343      if (FALSE && (average->objID == 0x1) && (average->catID)) {
     344        fprintf (stderr, "test obj\n");
     345      }
     346
     347      // mark the measurements matching this ranking
     348      markMeasureByRanking (&results->psfData[Nsec], measure, minRank, ID_MEAS_PHOTOM_PSF);
     349
     350      // when running -averages, we have no information about the images, so we cannot set this
     351      if (results->minUbercalDist[Nsec] > -1) {
     352        secfilt[Nsec].ubercalDist = results->minUbercalDist[Nsec];
     353      }
     354
     355      if (Nranking) {
     356        secfilt[Nsec].Mstdev = psfstats->sigma;
     357        secfilt[Nsec].Nused  = psfstats->Nmeas;
     358        secfilt[Nsec].Mmax   = psfstats->max;
     359        secfilt[Nsec].Mmin   = psfstats->min;
     360      }
     361
     362      secfilt[Nsec].psfQfMax     = results->psfQfMax[Nsec];
     363      secfilt[Nsec].psfQfPerfMax = results->psfQfPerfMax[Nsec];
     364
     365      if (isfinite(secfilt[Nsec].psfQfMax)     && (secfilt[Nsec].psfQfMax     > psfQfMax))     psfQfMax     = secfilt[Nsec].psfQfMax;
     366      if (isfinite(secfilt[Nsec].psfQfPerfMax) && (secfilt[Nsec].psfQfPerfMax > psfQfPerfMax)) psfQfPerfMax = secfilt[Nsec].psfQfPerfMax;
     367
     368      StatType *apstats = &results->apstats;
     369      Nranking = magStatsByRanking (&results->aperData[Nsec], apstats);
     370      if (Nranking) {
     371        secfilt[Nsec].Map     = apstats->mean;
     372        secfilt[Nsec].dMap    = apstats->error;
     373        secfilt[Nsec].sMap    = apstats->sigma;
     374        secfilt[Nsec].NusedAp = Nranking;
     375      }
     376      markMeasureByRanking (&results->aperData[Nsec], measure, minRank, ID_MEAS_PHOTOM_APER);
     377
     378      StatType *kronstats = &results->kronstats;
     379      Nranking = magStatsByRanking (&results->kronData[Nsec], kronstats);
     380      if (Nranking) {
     381        secfilt[Nsec].Mkron     = kronstats->mean;
     382        secfilt[Nsec].dMkron    = kronstats->error;
     383        secfilt[Nsec].sMkron    = kronstats->sigma;
     384        secfilt[Nsec].NusedKron = Nranking;
     385      }
     386      markMeasureByRanking (&results->kronData[Nsec], measure, minRank, ID_MEAS_PHOTOM_KRON);
     387
     388      // does this object appear extended in > 50% of measurements?
     389      if ((results->Next[Nsec] > 0) && (results->Next[Nsec] > 0.5*psfstats->Nmeas)) {
    394390        secfilt[Nsec].flags |= ID_SECF_OBJ_EXT;
    395391      }
    396392
    397       switch (pass) {
     393      switch (minRank) {
    398394        case 0:
    399           secfilt[Nsec].flags |= ID_PHOTOM_PASS_0;
    400           if (havePS1) GoodPS1 = TRUE;
     395          secfilt[Nsec].flags |= ID_SECF_RANK_0;
     396          if (results->havePS1[Nsec]) GoodPS1 = TRUE;
    401397          break;
    402398        case 1:
    403           secfilt[Nsec].flags |= ID_PHOTOM_PASS_1;
    404           if (havePS1) GoodPS1 = TRUE;
     399          secfilt[Nsec].flags |= ID_SECF_RANK_1;
     400          if (results->havePS1[Nsec]) GoodPS1 = TRUE;
    405401          break;
    406402        case 2:
    407           secfilt[Nsec].flags |= ID_PHOTOM_PASS_2;
    408           if (havePS1) GoodPS1 = TRUE;
     403          secfilt[Nsec].flags |= ID_SECF_RANK_2;
     404          if (results->havePS1[Nsec]) GoodPS1 = TRUE;
    409405          break;
    410406        case 3:
    411           secfilt[Nsec].flags |= ID_PHOTOM_PASS_3;
     407          secfilt[Nsec].flags |= ID_SECF_RANK_3;
    412408          break;
    413409        case 4:
    414           secfilt[Nsec].flags |= ID_PHOTOM_PASS_4;
     410          secfilt[Nsec].flags |= ID_SECF_RANK_4;
    415411          break;
    416       }
    417       if (haveSynth) {
    418         secfilt[Nsec].flags |= ID_SECF_USE_SYNTH;
    419       }
    420       if (havePS1) {
     412        default:
     413          break;
     414      }
     415      if (results->havePS1[Nsec]) {
    421416        secfilt[Nsec].flags |= ID_SECF_HAS_PS1;
    422417      }
    423       if (haveUbercal) {
     418      if (results->haveUbercal[Nsec]) {
    424419        secfilt[Nsec].flags |= ID_SECF_USE_UBERCAL;
    425       } 
     420      }
    426421    }
    427422  }
    428423
    429424  if (isSetMrelFinal) {
    430     if (pass == 0) {
    431       DVOAverageFlags flagBits = ID_OBJ_EXT | ID_OBJ_EXT_ALT | ID_OBJ_GOOD | ID_OBJ_GOOD_ALT;
    432       // we attempt to set a few flags here; reset those bits before trying:
    433       average[0].flags &= ~flagBits;
    434     }
    435 
    436425    if (NextPS1 && (NextPS1 > NpsfPS1)) {
    437426      average[0].flags |= ID_OBJ_EXT;
    … …  
    446435      average[0].flags |= ID_OBJ_GOOD_ALT;
    447436    }
     437    average[0].psfQF     = psfQfMax;
     438    average[0].psfQFperf = psfQfPerfMax;
    448439  }
    449440  return (TRUE);
    … …  
    455446// 2) select the BEST detections per filter (regardless of PRIMARY)
    456447// 3) apply the zero point and AB->Jy transformations
    457 int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t Nmeasure) {
     448int setMrelAverageStack (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr) {
     449
     450  // we are guaranteed to have average, measure, secfilt
     451  Average *average = &catalog[0].average[ave];
     452  off_t measureOffset = average->measureOffset;
     453  int Nmeasure = average->Nmeasure;
     454
     455  Measure *measure = &catalog[0].measure[measureOffset];
     456  SecFilt *secfilt = &catalog[0].secfilt[ave*Nsecfilt];
    458457
    459458  off_t k;
    460459
    461   float Mcal= 0, Mmos = 0, Mgrid = 0, Finst = 0;
     460  float Mcal = 0, Mmos = 0, Mgrid = 0, Finst = 0;
    462461
    463462  // set the primary projection cell and skycell for this coordinate
    … …  
    500499
    501500    off_t meas = measureOffset;
    502     for (k = 0; k < average[0].Nmeasure; k++, meas++) {
     501    for (k = 0; k < Nmeasure; k++, meas++) {
    503502
    504503      // only examine gpc1 stack data
    … …  
    528527        measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY;
    529528        secfilt[Nsec].stackPrmryOff = meas;
    530         myAssert (secfilt[Nsec].stackPrmryOff <= Nmeasure, "stackPrmryOff out of range");
     529        myAssert (secfilt[Nsec].stackPrmryOff <= catalog[0].Nmeasure, "stackPrmryOff out of range");
    531530      }
    532531
    … …  
    555554      }
    556555
    557       // NOTE: negative and insignificant vlues values are allowed, but not NAN flux values
     556      // NOTE: negative and insignificant flux values are allowed, but not NAN flux values
    558557      Finst = PhotFluxInst (&measure[k], MAG_CLASS_PSF);
    559558      if (isnan(Finst)) SKIP_THIS_MEAS_STACK(Ninst);
    560559
     560      // choose the best psfQFperf value for the BEST measurement
    561561      if (measure[k].psfQFperf > psfQFbest) {
    562562        psfQFbest = measure[k].psfQFperf;
     563        stackBestMeasure = k;
     564      }
     565      // ... UNLESS psfQFperf > 0.98 for the primary, in which case just use the primary.
     566      if ((measure[k].dbFlags & ID_MEAS_STACK_PRIMARY) && (measure[k].psfQFperf > 0.98)) {
    563567        stackBestMeasure = k;
    564568      }
    … …  
    640644
    641645    secfilt[Nsec].stackBestOff = k + measureOffset;
    642     myAssert (secfilt[Nsec].stackBestOff <= Nmeasure, "stackBestOff out of range");
     646    myAssert (secfilt[Nsec].stackBestOff <= catalog[0].Nmeasure, "stackBestOff out of range");
    643647
    644648    secfilt[Nsec].Nstack    = Nstack;
    … …  
    661665    average[0].flags |= ID_OBJ_BAD_STACK;
    662666  }
     667
    663668  return (TRUE);
    664669}
    665 
    666 # undef SKIP_THIS_MEAS
    667 # define SKIP_THIS_MEAS(REASON) {                               \
    668     measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM;  \
    669     results->REASON ++;                                         \
    670     continue; }
    671670
    672671// set mean of forced-warp measurements (selected by photcode range for now):
    673672// somewhat simplified relative to chip-photometry:
    674673// * no grid, no mosaic, no 2MASS, no SYNTH, no Ubercal, no flatcorr
    675 int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found, off_t Nmeasure) {
    676 
    677   off_t k;
    678   float Mcal= 0;
    679 
    680   // we are measuring means for 3 types of FLUXes: psf, ap, kron.  I am using the psf mag
    681   // error for the ap mags, but krons have their own errors.  it is an open question if I
    682   // should be doing weighted or unweighted fits (this is a user option)
    683   double *Fpsflist  = results->Mpsflist;
    684   double *dpsflist  = results->dpsflist;
    685   double *wpsflist  = results->wpsflist;
    686 
    687   double *Faplist   = results->Maplist;
    688   double *daplist   = results->daplist;
    689   double *waplist   = results->waplist;
    690 
    691   double *Fkronlist = results->Mkronlist;
    692   double *dkronlist = results->dkronlist;
    693   double *wkronlist = results->wkronlist;
    694 
    695   StatType *psfstats  = &results->psfstats;
    696   StatType *apstats   = &results->apstats;
    697   StatType *kronstats = &results->kronstats;
     674// analysis is done on flux, not mags (as the faintest objects will be nearly insignificant)
     675int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results) {
     676
     677  // we are guaranteed to have average, measure, secfilt
     678  Average *average    = &catalog[0].average[ave];
     679  off_t measureOffset = average->measureOffset;
     680
     681  Measure *measure    = &catalog[0].measure[measureOffset];
     682  SecFilt *secfilt    = &catalog[0].secfilt[ave*Nsecfilt];
     683  char *measureRank   = &catalog[0].measureRank[measureOffset];
     684
     685  off_t k, Nsec;
     686  float Mcal = 0;
     687
     688  SetMrelInfoResetObject (results); // reset the per-object arrays
    698689
    699690  // option for a test print
    700   if ((average[0].objID == 0xd6e) && (average[0].catID == 0x4984)) {
     691  if (FALSE && (average[0].objID == 0x0000020d) && (average[0].catID == 0x000049e1)) {
    701692    fprintf (stderr, "test obj\n");
    702693    print_measure_set_alt (average, secfilt, measure);
    703694  }
    704695
    705   // for a diff db, we just want primary, if available or not
    706   int priTessID, priProjID, priSkycellID;
     696  // for a diff db, we just want primary, if available or not. for warp we want to use the
     697  // skycell which provided the stack best entry
     698  int tessID, projID, skycellID;
    707699  if (IS_DIFF_DB) {
    708     get_tess_ids(&priTessID, &priProjID, &priSkycellID, average[0].R, average[0].D);
    709   }
    710 
    711   int Ns;
    712   for (Ns = 0; Ns < Nphotcodes; Ns++) {
    713 
    714     int thisCode = photcodes[Ns][0].code;
    715     int Nsec = GetPhotcodeNsec(thisCode);
    716 
    717     /* calculate the average mag in this SEC photcode for a single star */
    718 
    719     /* star/photcodes already calibrated */
    720     if (found[Nsec]) continue; 
    721      
    722     // I need to restrict the measurements coming only from the same skycell as the secfilt values
    723     // if we cannot make that association, skip this set of warps
    724     int stkTessID, stkProjID, stkSkycellID;
    725 
     700    get_tess_ids(&tessID, &projID, &skycellID, average[0].R, average[0].D);
     701  }
     702  for (Nsec = 0; Nsec < Nsecfilt; Nsec++) {
    726703    if (!IS_DIFF_DB) {
    727       myAssert (secfilt[Nsec].stackBestOff <= Nmeasure, "stackBestOff out of range");
    728       if (!FindImageSkycellID (secfilt[Nsec].stackBestOff, cat, &stkTessID, &stkProjID, &stkSkycellID)) continue;
    729     }
    730 
    731     float psfQfMax = 0.0;
    732     float psfQfPerfMax = 0.0;
    733 
    734     off_t meas = measureOffset;
    735 
    736     int Nap = 0;
    737     int Npsf = 0;
    738     int Nkron = 0;
    739 
    740     int Nwarp = 0;
    741     int NwarpGood = 0;
    742     for (k = 0; k < average[0].Nmeasure; k++, meas++) {
    743 
    744       // only examine gpc1 forced-warp data
    745       if (!isGPC1warp(measure[k].photcode)) continue;
    746 
    747       // skip measurements that do not match the current photcode
    748       PhotCode *code = GetPhotcodebyCode (measure[k].photcode);
    749       if (!code) continue;
    750       if (code->equiv != thisCode) { continue; }
    751 
    752       Nwarp ++;
    753       if (measure[k].psfQFperf > 0.85) NwarpGood ++;
    754 
    755       measure[k].dbFlags &= ~ID_MEAS_WARP_USED;
    756 
    757       if (IS_DIFF_DB) {
    758         if (MatchImageSkycellID (meas, cat, priTessID, priProjID, priSkycellID)) {
    759           measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY;
    760         } else {
    761           continue;
    762         }
    763       } else {
    764         if (!MatchImageSkycellID (meas, cat, stkTessID, stkProjID, stkSkycellID)) continue;
    765       }
    766  
    767       if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad);
    768 
    769       if (getImageEntry (meas, cat) < 0) {
    770         // measurements without an image are either external reference photometry or
    771         // data for which the associated image has not been loaded (probably because of
    772         // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
    773         Mcal = measure[k].Mcal; // check that this is zero for loaded REF value
    774       } else {
    775         // use getMcal not getMcal_alt?
    776         Mcal  = getMcal_alt (meas, cat, NULL, measure[k].Xccd, measure[k].Yccd);
    777         // Mcal  = getMcal (meas, cat);
    778         if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
    779       }
    780       float Fcal = MagToFlux(-Mcal);
    781 
    782       // these are only used below if IS_DIFF_DB
    783       if (IS_DIFF_DB) {
    784         if (measure[k].psfQF     > psfQfMax)     psfQfMax     = measure[k].psfQF;
    785         if (measure[k].psfQFperf > psfQfPerfMax) psfQfPerfMax = measure[k].psfQFperf;
    786       }
    787 
    788       // in the calculations below,
    789       // ...list gives the error per measurement, wlist gives the weight
    790       // we can modify the error and weight in a few ways:
    791       // 1) MIN_ERROR guarantees a floor
    792       // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
    793 
    794       // skip some absurd values NAN
    795       // NOTE : I am using PhotCat not PhotSys for now since GPC1 chip-to-chip color terms
    796       // are small (and not measured)
    797       float Fpsf = PhotFluxCat (&measure[k], MAG_CLASS_PSF);
    798       if (isnan(Fpsf)) SKIP_THIS_MEAS(Nsys);
    799       // if (Msys <  0.0) SKIP_THIS_MEAS(Nsys);
    800       // if (Msys > 30.0) SKIP_THIS_MEAS(Nsys);
    801 
    802       float dFpsf = PhotFluxCatErr (&measure[k], MAG_CLASS_PSF);
    803       if (isnan(dFpsf)) SKIP_THIS_MEAS (Nsys);
    804 
    805       int isBad = (measure[k].psfQF < 0.85);
    806       isBad |= isnan(measure[k].psfQF);
    807       if (isBad) SKIP_THIS_MEAS (Nsys);
    808 
    809       int isPoor = (measure[k].psfQFperf < 0.85);
    810       if ((pass == 0) && isPoor) SKIP_THIS_MEAS (Nsys);
    811 
    812       dFpsf = MAX (dFpsf, MIN_ERROR*Fpsf); // MIN_ERROR is a fractional error
    813       Fpsflist[Npsf] = Fpsf * Fcal;
    814       dpsflist[Npsf] = dFpsf * Fcal;
    815       wpsflist[Npsf] = 1.0;
    816       Npsf ++;
    817 
    818       measure[k].dbFlags |= ID_MEAS_WARP_USED;
    819 
    820       float Fap = PhotFluxCat (&measure[k], MAG_CLASS_APER);
    821       float dFap = PhotFluxCatErr (&measure[k], MAG_CLASS_APER);
    822       if (!isnan(Fap)) {
    823         Faplist[Nap] = Fap * Fcal;
    824         daplist[Nap] = dFap * Fcal;
    825         waplist[Nap] = 1.0; // drop weight lists?
    826         Nap ++;
    827       }
    828 
    829       float Fkron = PhotFluxCat (&measure[k], MAG_CLASS_KRON);
    830       float dFkron = PhotFluxCatErr (&measure[k], MAG_CLASS_KRON);
    831       if (!isnan(Fkron)) {
    832         Fkronlist[Nkron] = Fkron * Fcal;
    833         dkronlist[Nkron] = dFkron * Fcal;
    834         wkronlist[Nkron] = 1.0;
    835         Nkron ++;
    836       }
    837     }
    838     if (Npsf < 1) {
    839       secfilt[Nsec].Nwarp        = Nwarp;
    840       secfilt[Nsec].NwarpGood    = NwarpGood;
    841       secfilt[Nsec].NusedWrp     = Npsf;
    842       secfilt[Nsec].NusedApWrp   = Nap;
    843       secfilt[Nsec].NusedKronWrp = Nkron;
     704      myAssert (secfilt[Nsec].stackBestOff <= catalog->Nmeasure, "stackBestOff out of range");
     705      if (!FindImageSkycellID (secfilt[Nsec].stackBestOff, cat, &tessID, &projID, &skycellID)) continue;
     706    }
     707    results->tessID[Nsec]    = tessID;
     708    results->projID[Nsec]    = projID;
     709    results->skycellID[Nsec] = skycellID;
     710  }
     711
     712  // assign measurements to the photcode lists
     713  for (k = 0; k < average->Nmeasure; k++) {
     714    off_t meas = measureOffset + k;
     715
     716    measure[k].dbFlags &= ~ID_MEAS_WARP_USED; // this should not be set unless we use it in this pass
     717
     718    // skip measurements that do not have a valid photcode (raise exception?)
     719    PhotCode *code = GetPhotcodebyCode (measure[k].photcode);
     720    myAssert (code, "invalid photcode??");
     721
     722    // only examine gpc1 forced-warp data
     723    if (!isGPC1warp(measure[k].photcode)) continue;
     724
     725    // assign the Nsec value so we can assign to the right lists
     726    int Nsec = GetPhotcodeNsec (code->equiv);
     727    if (Nsec < 0) continue; // skip measurements which do not have an equiv average photcode
     728
     729    results->Nmeas[Nsec] ++;
     730    if (measure[k].psfQFperf > 0.85) results->NmeasGood[Nsec] ++;
     731
     732    // use primary skycell for DIFF, stack skycell for WARP
     733    if (MatchImageSkycellID (meas, cat, results->tessID[Nsec], results->projID[Nsec], results->skycellID[Nsec])) {
    844734      continue;
    845735    }
    846736
    847     found[Nsec] = TRUE;
    848 
    849     liststats (Fpsflist, dpsflist, wpsflist, Npsf, psfstats);
     737    // assign the max psfQF values for each Nsec value
     738    if (IS_DIFF_DB) {
     739      measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY;
     740      if (measure[k].psfQF     > results->psfQfMax[Nsec])     results->psfQfMax[Nsec]     = measure[k].psfQF;
     741      if (measure[k].psfQFperf > results->psfQfPerfMax[Nsec]) results->psfQfPerfMax[Nsec] = measure[k].psfQFperf;
     742    }
     743
     744    if (getImageEntry (meas, cat) < 0) {
     745      // measurements without an image are either external reference photometry or
     746      // data for which the associated image has not been loaded (probably because of
     747      // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
     748      Mcal = measure[k].Mcal; // check that this is zero for loaded REF value
     749    } else {
     750      // use getMcal not getMcal_alt?
     751      Mcal  = getMcal_alt (meas, cat, NULL, measure[k].Xccd, measure[k].Yccd);
     752      // Mcal  = getMcal (meas, cat);
     753      if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
     754    }
     755    float Fcal = MagToFlux(-Mcal);
     756
     757    // NOTE : I am using PhotCat not PhotSys for now since GPC1 chip-to-chip color terms
     758    // are small (and not measured)
     759    float Fpsf = PhotFluxCat (&measure[k], MAG_CLASS_PSF);
     760    float dFpsf = PhotFluxCatErr (&measure[k], MAG_CLASS_PSF);
     761    if (CHECK_VALID_FLUX(Fpsf, dFpsf)) {
     762      dFpsf = MAX (hypot(dFpsf, code->photomErrSys*Fpsf), MIN_ERROR*Fpsf); // bump up the error by a systematic floor
     763      int Npsf = results->psfData[Nsec].Nlist;
     764      results->psfData[Nsec].flxlist[Npsf] = Fpsf * Fcal;
     765      results->psfData[Nsec].errlist[Npsf] = dFpsf * Fcal;
     766      results->psfData[Nsec].wgtlist[Npsf] = 1.0;
     767      results->psfData[Nsec].ranking[Npsf] = measureRank ? measureRank[k] : 0;
     768      results->psfData[Nsec].measSeq[Npsf] = k;
     769      results->psfData[Nsec].Nlist ++;
     770    }     
     771
     772    float Fap = PhotFluxCat (&measure[k], MAG_CLASS_APER);
     773    float dFap = PhotFluxCatErr (&measure[k], MAG_CLASS_APER);
     774    if (CHECK_VALID_FLUX(Fap, dFap)) {
     775      dFap = MAX (hypot(dFap, code->photomErrSys*Fap), MIN_ERROR*Fap); // bump up the error by a systematic floor
     776      int Naper = results->aperData[Nsec].Nlist;
     777      results->aperData[Nsec].flxlist[Naper] = Fap * Fcal;
     778      results->aperData[Nsec].errlist[Naper] = dFap * Fcal;
     779      results->aperData[Nsec].wgtlist[Naper] = 1.0;
     780      results->aperData[Nsec].ranking[Naper] = measureRank ? measureRank[k] : 0;
     781      results->aperData[Nsec].measSeq[Naper] = k;
     782      results->aperData[Nsec].Nlist ++;
     783    }
     784
     785    float Fkron = PhotFluxCat (&measure[k], MAG_CLASS_KRON);
     786    float dFkron = PhotFluxCatErr (&measure[k], MAG_CLASS_KRON);
     787    if (CHECK_VALID_FLUX(Fkron, dFkron)) {
     788      dFkron = MAX (hypot(dFkron, code->photomErrSys*Fkron), MIN_ERROR*Fkron); // bump up the error by a systematic floor
     789      int Nkron = results->kronData[Nsec].Nlist;
     790      results->kronData[Nsec].flxlist[Nkron] = Fkron * Fcal;
     791      results->kronData[Nsec].errlist[Nkron] = dFkron * Fcal;
     792      results->kronData[Nsec].wgtlist[Nkron] = 1.0;
     793      results->kronData[Nsec].ranking[Nkron] = measureRank ? measureRank[k] : 0;
     794      results->kronData[Nsec].measSeq[Nkron] = k;
     795      results->kronData[Nsec].Nlist ++;
     796    }
     797  }
     798
     799  // find max values across filters
     800  float psfQfMax     = 0.0;
     801  float psfQfPerfMax = 0.0;
     802
     803  // now calculate the mean stats for the Nsec bands.
     804  for (Nsec = 0; Nsec < Nsecfilt; Nsec++) {
     805    dvo_secfilt_init (&secfilt[Nsec], SECFILT_RESET_WARP);
    850806
    851807    if (IS_DIFF_DB) {
    852808      // for non DIFF_DB, these are set in setMrelAverageExposure
    853       secfilt[Nsec].psfQfMax     = psfQfMax;
    854       secfilt[Nsec].psfQfPerfMax = psfQfPerfMax;
    855     }
    856 
    857     secfilt[Nsec].FpsfWrp  = psfstats->mean;
    858     secfilt[Nsec].dFpsfWrp = psfstats->error;
    859     secfilt[Nsec].sFpsfWrp = psfstats->sigma; // Mstdev is in millimags (not enough space for more precision)
    860     secfilt[Nsec].NusedWrp = psfstats->Nmeas;
    861     secfilt[Nsec].MpsfWrp  = isnan(secfilt[Nsec].FpsfWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FpsfWrp); // 8.9 since flux is in Jy
    862 
    863     // NOTE : use the modified weight for apmags as well as psf mags
    864     liststats (Faplist, daplist, waplist, Nap, apstats);
    865     secfilt[Nsec].FapWrp     = Nap > 0 ? apstats->mean : NAN;
    866     secfilt[Nsec].dFapWrp    = Nap > 0 ? apstats->error : NAN;
    867     secfilt[Nsec].sFapWrp    = Nap > 0 ? apstats->sigma : NAN;
    868     secfilt[Nsec].NusedApWrp = Nap;
    869     secfilt[Nsec].MapWrp  = isnan(secfilt[Nsec].FapWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FapWrp); // 8.9 since flux is in Jy
    870 
    871     liststats (Fkronlist, dkronlist, wkronlist, Nkron, kronstats);
    872     secfilt[Nsec].FkronWrp     = Nkron > 0 ? kronstats->mean  : NAN;
    873     secfilt[Nsec].dFkronWrp    = Nkron > 0 ? kronstats->error : NAN;
    874     secfilt[Nsec].sFkronWrp    = Nkron > 0 ? kronstats->sigma : NAN;
    875     secfilt[Nsec].NusedKronWrp = Nkron;
    876     secfilt[Nsec].MkronWrp  = isnan(secfilt[Nsec].FkronWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FkronWrp); // 8.9 since flux is in Jy
    877 
    878     secfilt[Nsec].Nwarp     = Nwarp;
    879     secfilt[Nsec].NwarpGood = NwarpGood;
     809      secfilt[Nsec].psfQfMax     = results->psfQfMax[Nsec];
     810      secfilt[Nsec].psfQfPerfMax = results->psfQfPerfMax[Nsec];
     811
     812      if (isfinite(secfilt[Nsec].psfQfMax)     && (secfilt[Nsec].psfQfMax     > psfQfMax))     psfQfMax     = secfilt[Nsec].psfQfMax;
     813      if (isfinite(secfilt[Nsec].psfQfPerfMax) && (secfilt[Nsec].psfQfPerfMax > psfQfPerfMax)) psfQfPerfMax = secfilt[Nsec].psfQfPerfMax;
     814    }
     815
     816    // if too few valid measurements meet the minimum criteria, go to the next entry
     817    StatType *psfstats = &results->psfstats;
     818    int Nranking = magStatsByRanking (&results->psfData[Nsec], psfstats);
     819    if (Nranking) {
     820      secfilt[Nsec].FpsfWrp  = psfstats->mean;
     821      secfilt[Nsec].dFpsfWrp = psfstats->error;
     822      secfilt[Nsec].sFpsfWrp = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
     823      secfilt[Nsec].NusedWrp = psfstats->Nmeas;
     824      secfilt[Nsec].MpsfWrp  = isnan(secfilt[Nsec].FpsfWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FpsfWrp); // 8.9 since flux is in Jy
     825    }
     826    int minRank = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10;
     827    markMeasureByRanking (&results->psfData[Nsec], measure, minRank, ID_MEAS_WARP_USED);
     828
     829    // if too few valid measurements meet the minimum criteria, go to the next entry
     830    StatType *apstats = &results->apstats;
     831    Nranking = magStatsByRanking (&results->aperData[Nsec], apstats);
     832    if (Nranking) {
     833      secfilt[Nsec].FapWrp     = apstats->mean;
     834      secfilt[Nsec].dFapWrp    = apstats->error;
     835      secfilt[Nsec].sFapWrp    = (apstats->Nmeas > 1) ? apstats->chisq : NAN;
     836      secfilt[Nsec].NusedApWrp = apstats->Nmeas;
     837      secfilt[Nsec].MapWrp     = isnan(secfilt[Nsec].FapWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FapWrp); // 8.9 since flux is in Jy
     838    }
     839
     840    // if too few valid measurements meet the minimum criteria, go to the next entry
     841    StatType *kronstats = &results->kronstats;
     842    Nranking = magStatsByRanking (&results->kronData[Nsec], kronstats);
     843    if (Nranking) {
     844      secfilt[Nsec].FkronWrp     = kronstats->mean;
     845      secfilt[Nsec].dFkronWrp    = kronstats->error;
     846      secfilt[Nsec].sFkronWrp    = (kronstats->Nmeas > 1) ? kronstats->chisq : NAN;
     847      secfilt[Nsec].NusedKronWrp = kronstats->Nmeas;
     848      secfilt[Nsec].MkronWrp     = isnan(secfilt[Nsec].FkronWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FkronWrp); // 8.9 since flux is in Jy
     849    }
     850
     851    secfilt[Nsec].Nwarp     = results->Nmeas[Nsec];
     852    secfilt[Nsec].NwarpGood = results->NmeasGood[Nsec];
    880853  }
    881854  return (TRUE);
    882855}
    883856
    884 int setGlobalObjStats (Average *average, Measure *measure) {
    885 
    886   off_t k;
    887 
    888   float psfQFmax     = 0.0;
    889   float psfQFperfmax = 0.0;
    890   float stargalmax   = 0.0;
    891 
    892   for (k = 0; k < average[0].Nmeasure; k++) {
    893 
    894     if (isfinite(measure[k].psfQF)) {
    895       psfQFmax     = MAX (psfQFmax,     measure[k].psfQF);
    896     }
    897     if (isfinite(measure[k].psfQFperf)) {
    898       psfQFperfmax = MAX (psfQFperfmax, measure[k].psfQFperf);
    899     }
    900 
    901     // only count psfQF > 0.85 here
    902     if (isfinite(measure[k].extNsigma) && (measure[k].psfQF > 0.85)) {
    903       stargalmax   = MAX (stargalmax,   measure[k].extNsigma);
    904     }
    905   }
    906   average[0].psfQF     = psfQFmax;
    907   average[0].psfQFperf = psfQFperfmax;
    908   average[0].stargal   = stargalmax;
    909 
    910   return (TRUE);
     857int magStatsByRanking (StatDataSet *dataset, StatType *stats) {
     858
     859  liststats_init (stats);
     860
     861  if (dataset->Nlist == 0) return 0;
     862
     863  // we have a list of measurements for this Nsec value, along with errors, weights, and their rank
     864  // first, sort by the rank (increasing)
     865  sort_StatDataSet (dataset);
     866
     867  // find the values with the same minimum rank:
     868  int Nranking = 0;
     869  int minRank = dataset->ranking[0]; // MIN (5, result->psfData[Nsec].ranking[0]);  -- only allow rank 5 or <
     870
     871  int i;
     872  for (i = 0; (i < dataset->Nlist) && (dataset->ranking[i] == minRank); i++, Nranking++);
     873
     874  liststats (dataset->flxlist, dataset->errlist, dataset->wgtlist, Nranking, stats);
     875  return (Nranking);
    911876}
     877
     878int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags flags) {
     879  int i;
     880
     881  for (i = 0; i < dataset->Nlist; i++) {
     882    int k = dataset->measSeq[i];
     883    if (dataset->ranking[i] > minrank) {
     884      measure[k].dbFlags &= ~flags;
     885    } else {
     886      measure[k].dbFlags |= flags;
     887    }
     888  }
     889  return TRUE;
     890}
     891
     892/* sort a coordinate pair (X,Y) and the associated index (S) */
     893void sort_StatDataSet (StatDataSet *dataset) {
     894 
     895# define SWAPFUNC(A,B){ double dtmp; int itmp; \
     896    dtmp = dataset->flxlist[A]; dataset->flxlist[A] = dataset->flxlist[B]; dataset->flxlist[B] = dtmp;     \
     897    dtmp = dataset->errlist[A]; dataset->errlist[A] = dataset->errlist[B]; dataset->errlist[B] = dtmp;     \
     898    dtmp = dataset->wgtlist[A]; dataset->wgtlist[A] = dataset->wgtlist[B]; dataset->wgtlist[B] = dtmp;     \
     899    itmp = dataset->ranking[A]; dataset->ranking[A] = dataset->ranking[B]; dataset->ranking[B] = itmp;     \
     900    itmp = dataset->measSeq[A]; dataset->measSeq[A] = dataset->measSeq[B]; dataset->measSeq[B] = itmp;     \
     901  }
     902# define COMPARE(A,B)(dataset->ranking[A] < dataset->ranking[B])
     903 
     904  OHANA_SORT (dataset->Nlist, COMPARE, SWAPFUNC);
     905 
     906# undef SWAPFUNC
     907# undef COMPARE
     908
     909}
     910

  • trunk/Ohana/src/relphot/src/setMrelFinal.c

    r37037 r38466  
    1212  int Nsecfilt = GetPhotcodeNsecfilt ();
    1313
    14   /* if we reset the catalog, reset all the current measurements */
     14  /*** RESET photometry flags and Mcal values ***/
    1515  if (RESET) {
    1616    // flags used by the photometry analysis (excluding UBERCAL)
    … …  
    4646      for (i = 0; i < catalog[0].Naverage; i++) {
    4747        off_t N = Nsecfilt*i+Nsec;
    48         dvo_secfilt_init (&catalog[0].secfilt[N]);
     48        dvo_secfilt_init (&catalog[0].secfilt[N], SECFILT_RESET_ALL);
    4949
    5050        off_t m = catalog[0].average[i].measureOffset;
    … …  
    100100  // XXX make this optional? (do not clean for -averages?)
    101101  if (!simpleAverage) clean_measures (catalog, 1, TRUE, flatcorr);    /* mark outliers ID_MEAS_POOR_PHOTOM */
    102   for (i = 0; i < 5; i++) {
    103     skip_measurements (catalog, i, flatcorr);       /* set ID_MEAS_SKIP for measures to be skipped */
    104     setMrelOutput  (catalog, 1, i, flatcorr);    /* set Mrel using remaining measures */
    105   }
     102
     103  ALLOCATE (catalog[0].measureRank, char, catalog[0].Nmeasure);
     104  setMeasureRank (catalog);
     105  setMrelOutput (catalog, 1, flatcorr);
    106106  if (!simpleAverage) setMcalOutput (catalog, 1, flatcorr);
    107107
    … …  
    267267}
    268268
    269 /* dophot types:
    270 
    271 1 - star
    272 2 - galaxy
    273 3 - star
    274 4 - weak
    275 5 - edge
    276 7 - weak
    277 9 - weak
    278 10 - sat
    279 */
    280 
     269void setMeasureRank (Catalog *catalog) {
     270
     271  int i;
     272
     273  // make these global
     274  int IMAGE_BAD  = ID_IMAGE_PHOTOM_NOCAL;
     275  int IMAGE_POOR = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
     276  int MEAS_BAD   = ID_MEAS_NOCAL        | ID_MEAS_SKIP_PHOTOM;
     277  int MEAS_POOR  = ID_MEAS_POOR_PHOTOM  | ID_MEAS_AREA;
     278
     279  Measure     *measure     = catalog[0].measure;
     280  MeasureTiny *measureT    = catalog[0].measureT;
     281  char        *measureRank = catalog[0].measureRank;
     282
     283  /* set measureRank[] based on various quality measurements */
     284  for (i = 0; i < catalog[0].Nmeasure; i++) {
     285    measureRank[i] = 11; // start at a low rank
     286
     287    // clear SKIP for all measures at first
     288    measureT[i].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
     289    measure [i].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
     290
     291    // skip measurements without a valid photcode
     292    PhotCode *code = GetPhotcodebyCode (measure[i].photcode);
     293    if (!code) continue;
     294
     295    // measurements outside time range have poor rank
     296    if (TimeSelect) {
     297      if (measure[i].t < TSTART)                     { measureRank[i] = 10; continue; }
     298      if (measure[i].t > TSTOP)                      { measureRank[i] = 10; continue; }
     299    }
     300   
     301    // measurements from BAD images and mosaics (not REF mags)
     302    off_t Nim = getImageEntry (i, 0);
     303    int imageFlags = getImageFlags (i, 0);
     304    int mosaicFlags = getMosaicFlags (i, 0);
     305
     306    if (Nim > -1) {
     307      if (imageFlags  & IMAGE_BAD)                   { measureRank[i] = 9; continue; }
     308
     309      if (mosaicFlags & IMAGE_BAD)                   { measureRank[i] = 9; continue; }
     310   
     311      // measurements ranked by inst mag limit (not REF mags)
     312      if (ImagSelect) {
     313        float mag = PhotInst (&measure[i], MAG_CLASS_PSF);
     314        if (mag < ImagMin)                           { measureRank[i] = 8; continue; }
     315        if (mag > ImagMax)                           { measureRank[i] = 8; continue; }
     316      }
     317    }
     318
     319    // RANK 1 : BAD photFlags
     320    if (measure[i].dbFlags & MEAS_BAD)               { measureRank[i] = 7; continue; }
     321
     322    // RANK 3 : BAD photFlags (eg, SAT, CR), internal outliers
     323    if (measure[i].photFlags & code->photomBadMask)  { measureRank[i] = 6; continue; }
     324
     325    if (Nim > -1) {
     326      if (imageFlags  & IMAGE_POOR)                  { measureRank[i] = 5; continue; }
     327      if (mosaicFlags & IMAGE_POOR)                  { measureRank[i] = 5; continue; }
     328    }
     329
     330    // RANK 2 : psfQF value
     331    if (!isfinite(measure[i].psfQF) || measure[i].psfQF < 0.85) { measureRank[i] = 4; continue; }
     332       
     333    // RANK 1 : POOR photFlags
     334    if (measure[i].dbFlags & MEAS_POOR)              { measureRank[i] = 3; continue; }
     335
     336    // RANK 1 : POOR photFlags
     337    if (measure[i].photFlags & code->photomPoorMask) { measureRank[i] = 2; continue; }
     338   
     339    // RANK 1 : psfQFperfect value
     340    if (!isfinite(measure[i].psfQFperf) || measure[i].psfQFperf < 0.85) { measureRank[i] = 1; continue; }
     341       
     342    // RANK 0 : perfect measurement:
     343    measureRank[i] = 0;
     344  }
     345}

  • trunk/Ohana/src/relphot/src/synthetic_zpts.c

    r38062 r38466  
    77static char *extname[5] = {"map_g", "map_r", "map_i", "map_z", "map_y"};
    88
    9 SynthZeroPoints *SynthZeroPointsLoad (char *filename) {
     9static SynthZeroPoints *zpts = NULL;
     10
     11int SynthZeroPointsLoad (char *filename) {
    1012
    1113  FILE *f = fopen (filename, "r");
    1214  if (!f) {
    1315    fprintf (stderr, "ERROR: cannot open file %s\n", filename);
    14     return NULL;
     16    return FALSE;
    1517  }
    1618
    … …  
    2325  // map_y
    2426
    25   SynthZeroPoints *zpts = NULL;
    2627  ALLOCATE (zpts, SynthZeroPoints, 1);
    2728 
    … …  
    2930    if (VERBOSE) fprintf (stderr, "can't read header\n");
    3031    fclose (f);
    31     return NULL;
     32    return FALSE;
    3233  }
    3334
    … …  
    4950  zpts->Ny = zpts->matrix[0].Naxis[1];
    5051
     52  return TRUE;
     53}
     54
     55SynthZeroPoints *SynthZeroPointsGet () {
    5156  return zpts;
    5257}
    53 
Note: See TracChangeset for help on using the changeset viewer.