IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38293


Ignore:
Timestamp:
May 21, 2015, 6:02:59 AM (11 years ago)
Author:
eugene
Message:

modifying relphot to better handle mean photomety assignment

Location:
branches/eam_branches/ipp-20150419/Ohana/src
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20150419/Ohana/src/libdvo/include/dvo.h

    r38154 r38293  
    117117  ID_MEAS_ICRF_QSO       = 0x00040000,  // this measurement is an ICRF reference position
    118118  ID_MEAS_IMAGE_EPOCH    = 0x00080000,  // this measurement is registered to the image epoch (not tied to ref catalog epoch)
     119  ID_MEAS_PHOTOM_PSF     = 0x00100000,  // this measurement is used for the mean psf mag
     120  ID_MEAS_PHOTOM_APER    = 0x00200000,  // this measurement is used for the mean ap mag
     121  ID_MEAS_PHOTOM_KRON    = 0x00400000,  // this measurement is used for the mean kron mag
    119122} DVOMeasureFlags;
    120123
    … …  
    786789  off_t *found_t;
    787790  off_t *foundWarp_t;
     791
     792  char *measureRank;
    788793//  off_t *image_t;
    789794//  off_t *mosaic_t;

  • branches/eam_branches/ipp-20150419/Ohana/src/relphot/include/relphot.h

    r38292 r38293  
    9999  double *errlist;            // mag errors for a star
    100100  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;
    101104} StatDataSet;
    102105
    … …  
    114117  // the final output steps.
    115118
     119  int Nsecfilt;
     120
    116121  StatDataSet  *psfData; // one is allocated for each primary (average) photcode
    117122  StatDataSet *aperData;
    … …  
    121126  double *psfqfperf_list;     // psfqfperf for all filters
    122127  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 *Next;                  // count of PS1 exposure (chip) measurements for this secfilt
     139  int *NexpPS1;               // count of PS1 exposure (chip) measurements for this secfilt
     140  int *haveUbercal;           // does this secfilt have any ubercal data?
     141
     142  float *minUbercalDist;
     143
     144  int NextPS1;                // count of PS1 exposure (chip) measurements for this secfilt
     145  int NpsfPS1;                // count of PS1 exposure (chip) measurements for this secfilt
    123146
    124147  StatType psfstats;
    … …  
    473496int print_measure_set_alt (Average *average, SecFilt *secfilt, Measure *measure);
    474497int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt);
    475 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);
     498int setMrelAverageExposure (off_t measureOffset, int cat, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT);
    476499int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t Nmeasure);
    477500int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found, off_t Nmeasure);
    … …  
    480503
    481504void SetMrelInfoInit (SetMrelInfo *results, int allocLists);
     505void SetMrelInfoFree (SetMrelInfo *results);
     506void SetMrelInfoReset (SetMrelInfo *results);
     507void SetMrelInfoResetObject (SetMrelInfo *results);
     508
     509void StatDataSetFree (StatDataSet *dataset, int Nsecfilt);
     510StatDataSet *StatDataSetAlloc (int Nsecfilt, int Nmax);
    482511
    483512int init_synthetic_mags ();
    … …  
    532561int isGPC1synth (int photcode);
    533562int whichGPC1filter (int photcode);
    534 
     563int is2MASS (int photcode);
    535564
    536565SynthZeroPoints *SynthZeroPointsLoad (char *filename);

  • branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/ImageOps.c

    r37807 r38293  
    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

  • branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/MosaicOps.c

    r37807 r38293  
    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;

  • branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/StarOps.c

    r38292 r38293  
    6464    ALLOCATE (dataset[i].wgtlist, double, Nmax);
    6565    ALLOCATE (dataset[i].errlist, double, Nmax);
     66    ALLOCATE (dataset[i].ranking,    int, Nmax);
     67    ALLOCATE (dataset[i].measSeq,    int, Nmax);
    6668  } 
    6769  return dataset;
    6870}
    6971
    70 void StatDataSetFree (StatDataSet *, int Nsecfilt) {
     72void StatDataSetFree (StatDataSet *dataset, int Nsecfilt) {
    7173
    7274  int i;
    … …  
    7678    FREE (dataset[i].wgtlist);
    7779    FREE (dataset[i].errlist);
     80    FREE (dataset[i].ranking);
     81    FREE (dataset[i].measSeq);
    7882  } 
    7983  FREE (dataset);
    … …  
    9094}
    9195
     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    results->minSYN[i]  = NAN;
     109
     110    results->psfQfMax[i]  = NAN;
     111    results->psfQfPerfMax[i]  = NAN;
     112
     113    results->Next[i] = 0;
     114    results->NexpPS1[i] = 0;
     115    results->havePS1[i] = 0;
     116    results->haveUbercal[i] = 0;
     117    results->minUbercalDist[i] = 1000;
     118
     119    results->psfData[i].Nlist = 0;
     120    results->aperData[i].Nlist = 0;
     121    results->kronData[i].Nlist = 0;
     122  }
     123  results->NextPS1 = 0;
     124  results->NpsfPS1 = 0;
     125}
     126
    92127void SetMrelInfoInit (SetMrelInfo *results, int allocLists) {
     128
     129  int i;
    93130  SetMrelInfoReset (results);
    94131  if (allocLists) {
    95132    results->Nsecfilt = GetPhotcodeNsecfilt ();
     133
    96134    results->psfData  = StatDataSetAlloc (results->Nsecfilt, Nmax);
    97135    results->aperData = StatDataSetAlloc (results->Nsecfilt, Nmax);
    98136    results->kronData = StatDataSetAlloc (results->Nsecfilt, Nmax);
     137
    99138    ALLOCATE (results->psfqf_list,     double, Nmax);
    100139    ALLOCATE (results->psfqfperf_list, double, Nmax);
    101140    ALLOCATE (results->stargal_list,   double, Nmax);
     141
     142    ALLOCATE (results->havePS1,   int, results->Nsecfilt);
     143    ALLOCATE (results->haveSYN,   int, results->Nsecfilt);
     144    ALLOCATE (results->needSYN,   int, results->Nsecfilt);
     145    ALLOCATE (results->measSYN,   int, results->Nsecfilt);
     146    ALLOCATE (results->minSYN,  float, results->Nsecfilt);
     147
     148    ALLOCATE (results->psfQfMax, float, results->Nsecfilt);
     149    ALLOCATE (results->psfQfPerfMax, float, results->Nsecfilt);
     150
     151    ALLOCATE (results->Next,           int, results->Nsecfilt);
     152    ALLOCATE (results->NexpPS1,        int, results->Nsecfilt);
     153    ALLOCATE (results->havePS1,        int, results->Nsecfilt);
     154    ALLOCATE (results->haveUbercal,    int, results->Nsecfilt);
     155
     156    ALLOCATE (results->minUbercalDist, float, results->Nsecfilt);
     157
     158    SetMrelInfoResetObject (results, results->Nsecfilt);
     159
     160    for (i = 0; i < NsynthNames; i++) {
     161      int photcode = GetPhotcodeCodebyName (synthNames[i]);
     162      if (!photcode) continue;
     163      int Nsec = GetPhotcodeNsec (photcode);
     164      if (Nsec < 0) continue;
     165      results->minSYN[Nsec] = synthMaxMags[i];
     166    }
    102167  } else {
    103168    results->Nsecfilt = 0;
    … …  
    105170    results->aperData = NULL;
    106171    results->kronData = NULL;
     172
    107173    results->psfqf_list     = NULL;
    108174    results->psfqfperf_list = NULL;
    109175    results->stargal_list   = NULL;
     176
     177    results->havePS1 = NULL;
     178    results->haveSYN = NULL;
     179    results->needSYN = NULL;
     180    results->measSYN = NULL;
     181    results->minSYN  = NULL;
     182
     183    results->psfQfMax = NULL;
     184    results->psfQfPerfMax = NULL;
     185    results->Next = NULL;
     186    results->NexpPS1 = NULL;
     187    results->havePS1 = NULL;
     188    results->haveUbercal = NULL;
     189    results->minUbercalDist = NULL;
    110190  }
    111191}
    112192
    113193void SetMrelInfoFree (SetMrelInfo *results) {
    114   StatDataSetFree (results->psfdata,  results->Nsecfilt);
    115   StatDataSetFree (results->aperdata, results->Nsecfilt);
    116   StatDataSetFree (results->krondata, results->Nsecfilt);
    117   FREE (results.psfqf_list);
    118   FREE (results.psfqfperf_list);
    119   FREE (results.stargal_list);
     194  StatDataSetFree (results->psfData,  results->Nsecfilt);
     195  StatDataSetFree (results->aperData, results->Nsecfilt);
     196  StatDataSetFree (results->kronData, results->Nsecfilt);
     197  FREE (results->psfqf_list);
     198  FREE (results->psfqfperf_list);
     199  FREE (results->stargal_list);
    120200}
    121201

  • branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/extra.c

    r38062 r38293  
    7070  return FALSE;
    7171}
     72
     73int is2MASS (int photcode) {
     74
     75  if ((photcode >= 2011) && (photcode <= 2013)) return TRUE;
     76  return FALSE;
     77}
     78

  • branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/liststats.c

    r34429 r38293  
    1414  fprintf (stderr, "ERROR: invalid stats mode: %s\n", strmode);
    1515  exit (1);
     16}
     17
     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;
    1631}
    1732

  • branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelCatalog.c

    r38292 r38293  
    22# define UBERCAL_WEIGHT 100.0
    33
    4 # define SKIP_THIS_MEAS(REASON) {                               \
    5     measureT[k].dbFlags |= ID_MEAS_SKIP_PHOTOM; \
     4# define SKIP_THIS_MEAS(REASON) {                       \
     5    measureT[k].dbFlags |= ID_MEAS_SKIP_PHOTOM;         \
    66    if (measure) {                                      \
    7       measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM;                \
    8     }                                                           \
    9     results->REASON ++;                                         \
     7      measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM;        \
     8    }                                                   \
     9    results->REASON ++;                                 \
    1010    continue; }
    1111
    … …  
    1313    measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM;          \
    1414    continue; }
     15
     16# define CHECK_VALID_MAG(MAG) (!isnan(MAG) && (MAG > 0.0) && (MAG < 30.0))
    1517
    1618static float MagToFlux (float Mag) {
    … …  
    5456    Measure *measure = catalog[Nc].measure ? &catalog[Nc].measure[m] : NULL;
    5557    SecFilt *secfilt = catalog[Nc].secfilt ? &catalog[Nc].secfilt[j*Nsecfilt] : NULL;
    56     setMrelAverageExposure (m, Nc, pass, flatcorr, results, average, &catalog[Nc].averageT[j], secfilt, measure, &catalog[Nc].measureT[m], &catalog[Nc].found_t[Nsecfilt*j]);
     58    setMrelAverageExposure (m, Nc, isSetMrelFinal, flatcorr, results, average, &catalog[Nc].averageT[j], secfilt, measure, &catalog[Nc].measureT[m], &catalog[Nc].found_t[Nsecfilt*j]);
    5759
    5860    // only apply Stack operation on setMrelFinal in first pass
    … …  
    6870    }
    6971  }
    70   // if (primaryCell) free (primaryCell);
     72
    7173  return (TRUE);
    7274}
    7375
     76// NOTE:
     77// Msys is measure[i].M + zp corrections
     78// Mcal is image[j].Mcal
     79// Mmos and Mgrid are offsets for mosaic and grid
     80
     81// tie down reference photometry if the -refcode (code) option is selected
     82// eg, -refcode g_SDSS
     83// this probably makes no sense in the context of multifilter analysis
     84// XXX probably need to use the photocde table to assign reference mag weights.
     85
     86// dlist gives the error per measurement, wlist gives the weight
     87// we can modify the error and weight in a few ways:
     88// 1) MIN_ERROR guarantees a floor
     89// 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
     90// 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
     91// 4) some reference photcode of some kind can be specified as fixed and have a high weight
     92
    7493// set mean of chip measurements (selected by photcode range for now):
    75 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) {
     94int setMrelAverageExposure (off_t measureOffset, int cat, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT) {
     95
     96  off_t k;
     97
     98  // we are measuring means for 3 types of mags: psf, ap, kron.
     99
     100  // in the final assignment, set the mean mags even if only 1 measurements exists
     101  int Nminmeas = isSetMrelFinal ? 1 : STAR_TOOFEW + 1;
     102
     103  // isSetMrelFinal : in the final pass, we set psf, kron, ap mag values and extra stats; other passes only do psf mags
     104
     105  // option for a test print
     106  if (FALSE && (averageT[0].objID == 0x7146) && (averageT[0].catID == 0x49d8)) {
     107    fprintf (stderr, "test obj\n");
     108    print_measure_set_alt (average, secfilt, measure);
     109  }
     110
     111  SetMrelInfoResetObject (results); // reset the per-object arrays
     112
     113  // assign measurements to the photcode lists
     114  for (k = 0; k < averageT[0].Nmeasure; k++) {
     115
     116    // skip measurements that do not have a valid photcode (raise exception?)
     117    PhotCode *code = GetPhotcodebyCode (measureT[k].photcode);
     118    myAssert (code, "invalid photcode??");
     119
     120    // SKIP gpc1 stack data (hard-wired photcodes)
     121    if (isGPC1stack(measureT[k].photcode)) continue;
     122   
     123    // SKIP gpc1 forced-warp data (hard-wired photcodes)
     124    if (isGPC1warp(measureT[k].photcode)) continue;
     125   
     126    if (isTycho(measureT[k].photcode)) { haveTYCHO = TRUE; }
     127
     128    if (is2MASS(measureT[k].photcode)) {
     129      if (measureT[k].photFlags & 0x00c00000) {
     130        Galaxy2MASS = TRUE; // per object value
     131      }
     132      if (measureT[k].photFlags & 0x00000007) {
     133        Good2MASS = TRUE; // per object value
     134      }
     135    }
     136
     137    // assign the Nsec value so we can assign to the right lists
     138    int Nsec = GetPhotcodeNsec (code->equiv);
     139    if (Nsec < 0) continue; // skip measurements which do not have an equiv average photcode
     140
     141    // force the use of SYN even if we have PS1 mags?
     142    if (isGPC1synth(measureT[k].photcode)) {
     143      results->haveSYN[Nsec] = TRUE;
     144      results->measSYN[Nsec] = k;
     145      if (measureT[k].M < results->minSYN[Nsec]) results->needSYN[Nsec] = TRUE;
     146      continue;
     147    }
     148
     149    // assign the max psfQF values for each Nsec value
     150    if (measureT[k].psfQF     > results->psfQfMax[Nsec])     results->psfQfMax[Nsec]     = measureT[k].psfQF;
     151    if (measureT[k].psfQFperf > results->psfQfPerfMax[Nsec]) results->psfQfPerfMax[Nsec] = measureT[k].psfQFperf;
     152   
     153    // are we a PS1 exposure photcode? (hard-wired photcodes)
     154    if (isGPC1chip(measureT[k].photcode)) { results->NexpPS1[Nsec] ++; results->havePS1[Nsec] = TRUE; }
     155
     156    // ** Choose the calibration (depends on the mode : do I have an image reference or not?)
     157    float Mcal = 0, Mmos = 0, Mgrid = 0;
     158    if (getImageEntry (measureOffset + k, cat) < 0) {
     159      // measurements without an image are either external reference photometry or
     160      // data for which the associated image has not been loaded (probably because of
     161      // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
     162      Mmos = Mgrid = 0;
     163      Mcal = measureT[k].Mcal; // check that this is zero for loaded REF value
     164    } else {
     165      Mcal  = getMcal_alt (meas, cat, flatcorr, measureT[k].Xccd, measureT[k].Yccd);
     166      if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
     167      Mmos  = getMmos  (meas, cat);
     168      if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
     169      Mgrid = getMgrid (meas, cat);
     170      if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
     171    }
     172
     173    // XXX ubercal dist should be by Nsec
     174    int myUbercalDist = getUbercalDist(meas, cat);
     175    result->minUbercalDist[Nsec] = MIN(result->minUbercalDist[Nsec], myUbercalDist);
     176
     177    int isUbercal = (measureT[k].dbFlags & ID_MEAS_PHOTOM_UBERCAL);
     178
     179    if (isUbercal) results->haveUbercal[Nsec] = TRUE; // haveUbercal is set per secfilt, isUbercal is per measure XXX define this array
     180
     181    int useUbercalWeight = isUbercal || (refPhotcode && (code->code == refPhotcode->code));
     182
     183    float Map, dMap;
     184    if (isSetMrelFinal) {
     185      Map = PhotCat (&measure[k], MAG_CLASS_APER);
     186      dMap = MAX (hypot(measureT[k].dMap,   code->photomErrSys), MIN_ERROR);
     187      if (CHECK_VALID_MAG(Map)) {
     188        // XXX need to carry a counter for each list set
     189        int Nap = results->apData[Nsec].Nlist;
     190        results->apData[Nsec].flxlist[Nap] = Map - Mcal - Mmos - Mgrid;
     191        results->apData[Nsec].errlist[Nap] = dMpsf; // no other option
     192        results->apData[Nsec].wgtlist[Nap] = useUbercalWeight ? UBERCAL_WEIGHT : 1.0;
     193        results->apData[Nsec].ranking[Nap] = measureRank[k];
     194        results->apData[Nsec].measSeq[Nap] = k;
     195        results->apData[Nsec].Nlist ++;
     196      }
     197    }
     198
     199    float Mkron, dMkron;
     200    if (isSetMrelFinal) {
     201      Mkron = PhotCat (&measure[k], MAG_CLASS_KRON);
     202      dMkron = MAX (hypot(measureT[k].dMkron, code->photomErrSys), MIN_ERROR);
     203      if (CHECK_VALID_MAG(Mkron)) {
     204        int Nkron = results->kronData[Nsec].Nlist;
     205        results->kronData[Nsec].flxlist[Nkron] = Mkron - Mcal - Mmos - Mgrid;
     206        results->kronData[Nsec].errlist[Nkron] = measure[k].dMkron;
     207        results->kronData[Nsec].wgtlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0;
     208        results->kronData[Nsec].ranking[Nkron] = measureRank[k];
     209        results->kronData[Nsec].measSeq[Nkron] = k;
     210        results->kronData[Nsec].Nlist ++;
     211      }
     212    }
     213
     214    float Mpsf = PhotSysTiny (&measureT[k], &averageT[0], &secfilt[0], MAG_CLASS_PSF);
     215    float dMpsf  = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR);
     216    if (CHECK_VALID_MAG(Mpsf)) {
     217      int Npsf = results->psfData[Nsec].Nlist;
     218      results->psfData[Nsec].flxlist[Npsf] = Mpsf - Mcal - Mmos - Mgrid;
     219      results->psfData[Nsec].errlist[Npsf] = dMpsf;
     220      results->psfData[Nsec].wgtlist[Npsf] = isUbercal ? UBERCAL_WEIGHT : 1.0;
     221      results->psfData[Nsec].ranking[Npsf] = measureRank[k];
     222      results->psfData[Nsec].measSeq[Npsf] = k;
     223      results->psfData[Nsec].Nlist ++;
     224    } else {
     225      measureT[k].dbFlags |= ID_MEAS_SKIP_PHOTOM; // XXX is this used outside of relphot?
     226      if (measure) measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM;
     227      results->Nsys ++;
     228    }
     229
     230    // count the extended detections (all PS1 bands)
     231    if (isGPC1chip(measureT[k].photcode) && !isnan(Map) && !isnan(Mpsf)) {
     232      float dMagAp = Mpsf - Map;
     233      float SigmaAp = hypot(0.1, 2.5*dMpsf);
     234      // XXX this is still quite ad hoc, but at least it:
     235      // (a) converges to 0.1 mag offset at the bright end
     236      // (b) converges to 0.5 mag offset at the faint end (dM = 0.2)
     237      if (dMagAp > SigmaAp) {
     238        results->Next[Nsec] ++;
     239        results->NextPS1 ++;
     240      } else {
     241        results->NpsfPS1 ++;
     242      }
     243    }
     244  }
     245
     246  // now calculate the mean stats for the Nsec bands.
     247  for (Nsec = 0; Nsec < Nsecfilt; Nsec++) {
     248    dvo_init_secfilt (&secfilt[Nsec]);
     249
     250    // XXX hardwired
     251    if ((thisCode < 6) || (thisCode == 9)) {
     252      secfilt[Nsec].Ncode = results->NexpPS1[Nsec];
     253    } else {
     254      secfilt[Nsec].Ncode = results->Ncode[Nsec]; // 2MASS data if it exists
     255    }
     256
     257    // XXX need to also supply ZP from map, apply tycho logic
     258    if (results->needSYN[Nsec]) {
     259      if (!isSetMrelFinal) continue;
     260
     261      // use the single SYNTH value instead of the other mags here
     262      float Mpsf = measure[results->measSYN[Nsec]].M;
     263      float ZP = 0.0;
     264
     265      if (!haveTycho) {
     266        // need to look up the (X,Y) coords from this (R,D) location
     267        double X, Y;
     268        double R = ohana_normalize_angle_to_midpoint (measure->R, 0.0); // XXX 0.0 or 180.0 for center?
     269        RD_to_XY (&X, &Y, R, measure->D, &zpts->coords);
     270        if (X < 0) return FALSE;
     271        if (Y < 0) return FALSE;
     272        if (X >= zpts->Nx) return FALSE;
     273        if (Y >= zpts->Ny) return FALSE;
     274       
     275        int Xpix = X;
     276        int Ypix = Y;
     277        int Npix = Xpix + Ypix*zpts->Nx;
     278       
     279        float *value = (float *) zpts->matrix[Nsec].buffer;
     280        ZP = !isnan(value[Npix]) ? value[Npix] : 0.0;
     281        secfilt[Nsec].flags |= ID_SECF_HAVE_TYCHO;
     282      }
     283      secfilt[Nsec].M      = Mpsf + ZP;
     284      secfilt[Nsec].dM     = 0.6;
     285      secfilt[Nsec].chipsq = 0.0;
     286      secfilt[Nsec].flags |= ID_SECF_USE_SYNTH;
     287      continue;
     288    }
     289
     290    // if too few valid measurements meet the minimum criteria, go to the next entry
     291    int Nranking = magStatsByRanking (&results->psfData[Nsec], &results->psfstats);
     292    if (Nranking < Nminmeas) {
     293      secfilt[Nsec].flags |= ID_STAR_FEW;
     294    } else {
     295      secfilt[Nsec].M    = psfstats->mean;
     296      secfilt[Nsec].dM   = psfstats->error;
     297      secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
     298    }
     299    int minRank = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10;
     300
     301    // when running -averages, we have no information about the images, so we cannot set this
     302    if (minUbercalDist > -1) {
     303      secfilt[Nsec].ubercalDist = minUbercalDist;
     304    }
     305
     306    if (isSetMrelFinal) {
     307      secfilt[Nsec].Mstdev = psfstats->sigma;
     308      secfilt[Nsec].Nused  = psfstats->Nmeas;
     309      secfilt[Nsec].Mmax = psfstats->max;
     310      secfilt[Nsec].Mmin = psfstats->min;
     311
     312      secfilt[Nsec].psfQfMax     = psfQfMax;
     313      secfilt[Nsec].psfQfPerfMax = psfQfPerfMax;
     314
     315      Nranking = magStatsByRanking (&results->aperData[Nsec], results->apstats);
     316      secfilt[Nsec].Map     = apstats->mean;
     317      secfilt[Nsec].dMap    = apstats->error;
     318      secfilt[Nsec].sMap    = apstats->sigma;
     319      secfilt[Nsec].NusedAp = Nranking;
     320
     321      Nranking = magStatsByRanking (&results->kronData[Nsec], kronstats);
     322      secfilt[Nsec].Mkron     = kronstats->mean;
     323      secfilt[Nsec].dMkron    = kronstats->error;
     324      secfilt[Nsec].sMkron    = kronstats->sigma;
     325      secfilt[Nsec].NusedKron = Nranking;
     326
     327      // does this object appear extended in > 50% of measurements?
     328      if ((results->Next[Nsec] > 0) && (results->Next[Nsec] > 0.5*results->Npsf[Nsec])) {
     329        secfilt[Nsec].flags |= ID_SECF_OBJ_EXT;
     330      }
     331
     332      switch (minRank) {
     333        case 0:
     334          secfilt[Nsec].flags |= ID_PHOTOM_PASS_0;
     335          if (results->havePS1[Nsec]) GoodPS1 = TRUE;
     336          break;
     337        case 1:
     338          secfilt[Nsec].flags |= ID_PHOTOM_PASS_1;
     339          if (results->havePS1[Nsec]) GoodPS1 = TRUE;
     340          break;
     341        case 2:
     342          secfilt[Nsec].flags |= ID_PHOTOM_PASS_2;
     343          if (results->havePS1[Nsec]) GoodPS1 = TRUE;
     344          break;
     345        case 3:
     346          secfilt[Nsec].flags |= ID_PHOTOM_PASS_3;
     347          break;
     348        case 4:
     349          secfilt[Nsec].flags |= ID_PHOTOM_PASS_4;
     350          break;
     351        default:
     352          break;
     353      }
     354      if (results->havePS1[Nsec]) {
     355        secfilt[Nsec].flags |= ID_SECF_HAS_PS1;
     356      }
     357      if (results->haveUbercal[Nsec]) {
     358        secfilt[Nsec].flags |= ID_SECF_USE_UBERCAL;
     359      }
     360    }
     361  }
     362
     363  if (isSetMrelFinal) {
     364    if (results->NextPS1 && (results->NextPS1 > results->NpsfPS1)) {
     365      average[0].flags |= ID_OBJ_EXT;
     366    }
     367    if (GoodPS1) {
     368      average[0].flags |= ID_OBJ_GOOD;
     369    }
     370    if (Galaxy2MASS) {
     371      average[0].flags |= ID_OBJ_EXT_ALT;
     372    }
     373    if (Good2MASS) {
     374      average[0].flags |= ID_OBJ_GOOD_ALT;
     375    }
     376  }
     377  return (TRUE);
     378}
     379
     380// set mean of chip measurements (selected by photcode range for now):
     381int setMrelAverageExposure_old (off_t measureOffset, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found) {
    76382
    77383  off_t k;
    … …  
    81387  // error for the ap mags, but krons have their own errors.  it is an open question if I
    82388  // should be doing weighted or unweighted fits (this is a user option)
    83 
    84   // these lists are pre-allocated to the max length
    85   double **Mpsflist  = results->Mpsflist;
    86   double **dpsflist  = results->dpsflist;
    87   double **wpsflist  = results->wpsflist;
    88 
    89   double **Maplist   = results->Maplist;
    90   double **daplist   = results->daplist;
    91   double **waplist   = results->waplist;
    92 
    93   double **Mkronlist = results->Mkronlist;
    94   double **dkronlist = results->dkronlist;
    95   double **wkronlist = results->wkronlist;
     389  double *Mpsflist  = results->Mpsflist;
     390  double *dpsflist  = results->dpsflist;
     391  double *wpsflist  = results->wpsflist;
     392
     393  double *Maplist   = results->Maplist;
     394  double *daplist   = results->daplist;
     395  double *waplist   = results->waplist;
     396
     397  double *Mkronlist = results->Mkronlist;
     398  double *dkronlist = results->dkronlist;
     399  double *wkronlist = results->wkronlist;
    96400
    97401  StatType *psfstats  = &results->psfstats;
    … …  
    105409    fprintf (stderr, "test obj\n");
    106410    print_measure_set_alt (average, secfilt, measure);
    107   }
    108 
    109   // first pass: check for PS1[n], Synth, Tycho, 2MASS
    110 
    111   // first, determine the state of things for each of the main filters:
    112   // in this bit, we are using hardwired photcode values:
    113   // PS1_g, PS1_r, PS1_i, PS1_z, PS1_y, PS1_w == (1,2,3,4,5,6); 0 == none
    114   // note that Nsec = PS_N - 1
    115 
    116   int havePS1[7];
    117   int haveSYN[7];
    118   int needSYN[7];
    119   int measSYN[7];
    120 
    121   int haveTYCHO = FALSE;
    122   int have2MASS = FALSE;
    123 
    124   // was: 13.64, 13.76, 13.74, 12.94, 12.01
    125   float MaxMagForceSynth[] = {0.0, 13.5, 13.5, 13.5, 13.0, 12.0, 13.5};
    126   //                          NONE, g,    r,    i,    z,    y,    w
    127 
    128   for (i = 0; i < 7; i++) {
    129     havePS1[i] = FALSE; // do we have any PS1 measurement for this filter?
    130     haveSYN[i] = FALSE; // do we have any synthetic measurement for this filter?
    131     needSYN[i] = FALSE; // force the use of the synthetic measurement for this filter?
    132     measSYN[i] = -1; // which entries carries the synthetic measurement for this filter?
    133   }
    134 
    135   for (k = 0; k < average[0].Nmeasure; k++) {
    136 
    137     // skip measurements that do not have a valid photcode (which ones lack photcode?)
    138     PhotCode *code = GetPhotcodebyCode (measure[k].photcode);
    139     if (!code) continue;
    140 
    141     // is this synth_grizy?
    142     // is this PS1?
    143 
    144     // are we a PS1 exposure photcode?
    145     if (isGPC1chip(measure[k].photcode)) {
    146       int Nfilter = whichGPC1filter (measure[k].photcode);
    147       havePS1[Nfilter] = TRUE;
    148     }
    149     if (isGPC1synth(measure[k].photcode)) {
    150       int Nfilter = measure[k].photcode - 3000;
    151       haveSYN[Nfilter] = TRUE;
    152       measSYN[Nfilter] = k;
    153 
    154       // force the use of SYN even if we have PS1 mags?
    155       if (measure[k].M < MaxMagForceSynth[Nfilter]) {
    156         needSYN[Nfilter] = TRUE;
    157       }
    158     }
    159     if (isTycho(measure[k].photcode)) {
    160       haveTYCHO = TRUE;
    161     }
    162     if (is2MASS(measure[k].photcode)) {
    163       have2MASS = TRUE;
    164     }
    165411  }
    166412
    … …  
    505751}
    506752
    507 // set mean of chip measurements (selected by photcode range for now):
    508 int setMrelAverageExposure_old (off_t measureOffset, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found) {
    509 
    510   off_t k;
    511   float Msys = 0, Mcal= 0, Mmos = 0, Mgrid = 0;
    512 
    513   // we are measuring means for 3 types of mags: psf, ap, kron.  I am using the psf mag
    514   // error for the ap mags, but krons have their own errors.  it is an open question if I
    515   // should be doing weighted or unweighted fits (this is a user option)
    516   double *Mpsflist  = results->Mpsflist;
    517   double *dpsflist  = results->dpsflist;
    518   double *wpsflist  = results->wpsflist;
    519 
    520   double *Maplist   = results->Maplist;
    521   double *daplist   = results->daplist;
    522   double *waplist   = results->waplist;
    523 
    524   double *Mkronlist = results->Mkronlist;
    525   double *dkronlist = results->dkronlist;
    526   double *wkronlist = results->wkronlist;
    527 
    528   StatType *psfstats  = &results->psfstats;
    529   StatType *apstats   = &results->apstats;
    530   StatType *kronstats = &results->kronstats;
    531 
    532   int isSetMrelFinal = (pass >= 0);
    533 
    534   // option for a test print
    535   if (FALSE && (average[0].objID == 0x7146) && (average[0].catID == 0x49d8)) {
    536     fprintf (stderr, "test obj\n");
    537     print_measure_set_alt (average, secfilt, measure);
    538   }
    539 
    540   int GoodPS1 = FALSE;
    541   int Good2MASS = FALSE;
    542   int Galaxy2MASS = FALSE;
    543 
    544   int NextPS1 = 0;
    545   int NpsfPS1 = 0;
    546 
    547   int Ns;
    548   for (Ns = 0; Ns < Nphotcodes; Ns++) {
    549 
    550     int thisCode = photcodes[Ns][0].code;
    551     int Nsec = GetPhotcodeNsec(thisCode);
    552 
    553     /* calculate the average mag in this SEC photcode for a single star */
    554 
    555     /* star/photcodes already calibrated */
    556     if ( isSetMrelFinal && found[Nsec]) continue; 
    557      
    558     // skip bad stars
    559     if (!isSetMrelFinal && (secfilt[Nsec].flags & STAR_BAD)) continue;
    560 
    561     int NexpPS1 = 0;
    562     int Ncode = 0;
    563     int Next = 0;
    564     int haveSynth = FALSE;
    565     int havePS1   = FALSE;
    566 
    567     int forceSynth = FALSE;
    568     int forceSynthEntry = -1;
    569 
    570     int haveUbercal = FALSE;
    571 
    572     int minUbercalDist = 1000;
    573    
    574     off_t meas = measureOffset;
    575 
    576     float psfQfMax = 0.0;
    577     float psfQfPerfMax = 0.0;
    578 
    579     int Nap = 0;
    580     int Npsf = 0;
    581     int Nkron = 0;
    582     for (k = 0; k < averageT[0].Nmeasure; k++, meas++) {
    583 
    584       // skip measurements that do not match the current photcode
    585       PhotCode *code = GetPhotcodebyCode (measureT[k].photcode);
    586       if (!code) continue;
    587       if (code->equiv != thisCode) { continue; }
    588       Ncode ++;
    589 
    590       // are we a PS1 exposure photcode?
    591       if (isGPC1chip(measureT[k].photcode)) NexpPS1 ++;
    592 
    593       // SKIP gpc1 stack data
    594       if (isGPC1stack(measureT[k].photcode)) continue;
    595 
    596       // SKIP gpc1 forced-warp data
    597       if (isGPC1warp(measureT[k].photcode)) continue;
    598 
    599       if (isSetMrelFinal) {
    600         if (measure[k].psfQF     > psfQfMax)     psfQfMax     = measure[k].psfQF;
    601         if (measure[k].psfQFperf > psfQfPerfMax) psfQfPerfMax = measure[k].psfQFperf;
    602       }
    603 
    604       if (measureT[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad);
    605 
    606       if (getImageEntry (meas, cat) < 0) {
    607         // measurements without an image are either external reference photometry or
    608         // data for which the associated image has not been loaded (probably because of
    609         // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
    610         Mmos = Mgrid = 0;
    611         Mcal = measureT[k].Mcal; // check that this is zero for loaded REF value
    612       } else {
    613         Mcal  = getMcal_alt (meas, cat, flatcorr, measureT[k].Xccd, measureT[k].Yccd);
    614         if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
    615         Mmos  = getMmos  (meas, cat);
    616         if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
    617         Mgrid = getMgrid (meas, cat);
    618         if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
    619       }
    620 
    621       // skip some absurd values NAN, < 0.0, > 30.0
    622       Msys = PhotSysTiny (&measureT[k], &averageT[0], &secfilt[0], MAG_CLASS_PSF);
    623       if (isnan(Msys)) SKIP_THIS_MEAS(Nsys);
    624       if (Msys <  0.0) SKIP_THIS_MEAS(Nsys);
    625       if (Msys > 30.0) SKIP_THIS_MEAS(Nsys);
    626 
    627       int myUbercalDist = getUbercalDist(meas, cat);
    628       minUbercalDist = MIN(minUbercalDist, myUbercalDist);
    629 
    630       int isUbercal = (measureT[k].dbFlags & ID_MEAS_PHOTOM_UBERCAL);
    631       float dMpsf = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR);
    632 
    633       if (isUbercal) haveUbercal = TRUE; // haveUbercal is set per secfilt, isUbercal is per measure
    634 
    635       if (isSetMrelFinal) {
    636         // special options for PS1 data
    637         if ((measure[k].photcode >= 10000) && (measure[k].photcode <= 10500)) {
    638           // count the extended detections
    639           if (!isnan(measure[k].Map)) {
    640             float dMagAp = measure[k].M - measure[k].Map;
    641             float SigmaAp = hypot(0.1, 2.5*measure[k].dM);
    642             // XXX this is still quite ad hoc, but at least it:
    643             // (a) converges to 0.1 mag offset at the bright end
    644             // (b) converges to 0.5 mag offset at the faint end (dM = 0.2)
    645             if (dMagAp > SigmaAp) {
    646               Next ++;
    647               NextPS1 ++;
    648             } else {
    649               NpsfPS1 ++;
    650             }
    651           }
    652           havePS1 = TRUE;
    653         }
    654 
    655         // count extended detections for 2MASS (XXX NOTE hardwired photcodes 2011, 2012, 2013)
    656         if ((measure[k].photcode >= 2011) && (measure[k].photcode <= 2013)) {
    657           if (measure[k].photFlags & 0x00c00000) {
    658             Next ++;
    659             Galaxy2MASS = TRUE;
    660           }
    661           if (pass == 0) {
    662             if (measure[k].photFlags & 0x00000007) {
    663               Good2MASS = TRUE;
    664             } else {
    665               // detections without one of these bits should only be used in PASS_1
    666               SKIP_THIS_MEAS(Nbad);
    667             }
    668           }
    669         }
    670 
    671         // Blindly accepth the SYNTH mags if we are above saturation, otherwise,
    672         // ignore SYNTH photcodes until PASS == 4 (where we also accept saturated stars)
    673         if ((measure[k].photcode >= 3001) && (measure[k].photcode <= 3005)) {
    674           // something of a hack: force object to use synth values if synth mags >>
    675           // saturation (3pi instrumental mags < -15)
    676           float MaxMagForceSynth = NAN;
    677           switch (measure[k].photcode) {
    678             case 3001:
    679               MaxMagForceSynth = 13.64;
    680               break;
    681             case 3002:
    682               MaxMagForceSynth = 13.76;
    683               break;
    684             case 3003:
    685               MaxMagForceSynth = 13.74;
    686               break;
    687             case 3004:
    688               MaxMagForceSynth = 12.94;
    689               break;
    690             case 3005:
    691               MaxMagForceSynth = 12.01;
    692               break;
    693           }
    694           if (measureT[k].M < MaxMagForceSynth) {
    695             forceSynth = TRUE;
    696             forceSynthEntry = Npsf;
    697           } else {
    698             if (pass < 4) {
    699               SKIP_THIS_MEAS(Nbad);
    700             }
    701             haveSynth = TRUE;
    702           }
    703         }
    704 
    705         // Map (Maplist) and Mkron (Mkronlist,dkronlist) are used to calculate mean mags per filter
    706         float Map = PhotCat (&measure[k], MAG_CLASS_APER);
    707         if (!isnan(Map)) {
    708             Maplist[Nap] = Map - Mcal - Mmos - Mgrid;
    709             daplist[Nap] = dMpsf;// XXX check on this...
    710             waplist[Nap] = isUbercal ? UBERCAL_WEIGHT : 1.0;
    711             Nap ++;
    712         }
    713 
    714         float Mkron = PhotCat (&measure[k], MAG_CLASS_KRON);
    715         if (!isnan(Mkron)) {
    716           Mkronlist[Nkron] = Mkron - Mcal - Mmos - Mgrid;
    717           dkronlist[Nkron] = measure[k].dMkron;
    718           wkronlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0;
    719           Nkron ++;
    720         }
    721       } // if (isSetMrelFinal)
    722 
    723       // dlist gives the error per measurement, wlist gives the weight
    724       // we can modify the error and weight in a few ways:
    725       // 1) MIN_ERROR guarantees a floor
    726       // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
    727       // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
    728       // 4) some reference photcode of some kind can be specified as fixed and have a high weight
    729       Mpsflist[Npsf] = Msys - Mcal - Mmos - Mgrid;
    730       dpsflist[Npsf] = dMpsf;
    731       wpsflist[Npsf] = isUbercal ? UBERCAL_WEIGHT : 1.0;
    732 
    733       // NOTE:
    734       // Msys is measure[i].M + zp corrections
    735       // Mcal is image[j].Mcal
    736       // Mmos and Mgrid are offsets for mosaic and grid
    737 
    738       // tie down reference photometry if the -refcode (code) option is selected
    739       // eg, -refcode g_SDSS
    740       // this probably makes no sense in the context of multifilter analysis
    741       // XXX probably need to use the photocde table to assign reference mag weights.
    742       if (refPhotcode) {
    743         if (code->code == refPhotcode->code) {
    744           wpsflist[Npsf] = UBERCAL_WEIGHT;
    745         }
    746       }
    747       Npsf ++;
    748     }
    749 
    750     int Nminmeas = isSetMrelFinal ? 1 : STAR_TOOFEW + 1;
    751 
    752     // XXX : ugh : another hard-wired photcode entry...
    753     if (isSetMrelFinal && (pass == 0)) {
    754       if ((thisCode < 6) || (thisCode == 9)) {
    755         secfilt[Nsec].Ncode = NexpPS1;
    756       } else {
    757         secfilt[Nsec].Ncode = Ncode; // 2MASS data if it exists
    758       }
    759     }
    760 
    761     // when performing the grid analysis, STAR_TOOFEW should be set to 1;
    762     if (Npsf < Nminmeas) { /* too few measurements */
    763       // fprintf (f, "%10.6f %10.6f %d %d %d\n", averageT[0].R, averageT[0].D, measureT[0].imageID, Npsf, STAR_TOOFEW);
    764       secfilt[Nsec].flags |= ID_STAR_FEW;
    765       if (Ncode == 0) {
    766         results->Ncode ++;
    767       } else {
    768         results->Nfew ++;
    769       }
    770       continue;
    771     } else {
    772       secfilt[Nsec].flags &= ~ID_STAR_FEW;
    773     }   
    774 
    775     if (forceSynth) {
    776       // use the single SYNTH value instead of the other mags here
    777       myAssert ((forceSynthEntry < Npsf) && (forceSynthEntry >= 0), "programming error");
    778       Mpsflist[0] = Mpsflist[forceSynthEntry];
    779       dpsflist[0] = dpsflist[forceSynthEntry];
    780       wpsflist[0] = wpsflist[forceSynthEntry];
    781       Npsf = 1;
    782     }
    783     liststats (Mpsflist, dpsflist, wpsflist, Npsf, psfstats);
    784 
    785     secfilt[Nsec].M      = psfstats->mean;
    786     secfilt[Nsec].dM     = psfstats->error;
    787     secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
    788 
    789     // when running -averages, we have no information about the images, so we cannot set this
    790     if (minUbercalDist > -1) {
    791       secfilt[Nsec].ubercalDist = minUbercalDist;
    792     }
    793 
    794     if (isSetMrelFinal) {
    795       found[Nsec] = TRUE;
    796 
    797       secfilt[Nsec].Mstdev = psfstats->sigma; // Mstdev is in millimags (not enough space for more precision)
    798       // secfilt[Nsec].Ncode = Ncode;
    799       secfilt[Nsec].Nused = psfstats->Nmeas;
    800 
    801       secfilt[Nsec].Mmax = psfstats->max;
    802       secfilt[Nsec].Mmin = psfstats->min;
    803 
    804       secfilt[Nsec].psfQfMax     = psfQfMax;
    805       secfilt[Nsec].psfQfPerfMax = psfQfPerfMax;
    806 
    807       // NOTE : use the modified weight for apmags as well as psf mags
    808       liststats (Maplist, daplist, waplist, Nap, apstats);
    809       secfilt[Nsec].Map  = Nap > 0 ? apstats->mean : NAN;
    810       secfilt[Nsec].dMap  = Nap > 0 ? apstats->error : NAN;
    811       secfilt[Nsec].sMap  = Nap > 0 ? apstats->sigma : NAN;
    812       secfilt[Nsec].NusedAp  = Nap;
    813 
    814       liststats (Mkronlist, dkronlist, wkronlist, Nkron, kronstats);
    815       secfilt[Nsec].Mkron  = Nkron > 0 ? kronstats->mean  : NAN;
    816       secfilt[Nsec].dMkron = Nkron > 0 ? kronstats->error : NAN;
    817       secfilt[Nsec].sMkron = Nkron > 0 ? kronstats->sigma : NAN;
    818       secfilt[Nsec].NusedKron  = Nkron;
    819 
    820       // NOTE: for 2MASS measurements, Next should be 1, as should N
    821       if ((Next > 0) && (Next > 0.5*Npsf)) {
    822         secfilt[Nsec].flags |= ID_SECF_OBJ_EXT;
    823       }
    824 
    825       switch (pass) {
    826         case 0:
    827           secfilt[Nsec].flags |= ID_PHOTOM_PASS_0;
    828           if (havePS1) GoodPS1 = TRUE;
    829           break;
    830         case 1:
    831           secfilt[Nsec].flags |= ID_PHOTOM_PASS_1;
    832           if (havePS1) GoodPS1 = TRUE;
    833           break;
    834         case 2:
    835           secfilt[Nsec].flags |= ID_PHOTOM_PASS_2;
    836           if (havePS1) GoodPS1 = TRUE;
    837           break;
    838         case 3:
    839           secfilt[Nsec].flags |= ID_PHOTOM_PASS_3;
    840           break;
    841         case 4:
    842           secfilt[Nsec].flags |= ID_PHOTOM_PASS_4;
    843           break;
    844       }
    845       if (haveSynth) {
    846         secfilt[Nsec].flags |= ID_SECF_USE_SYNTH;
    847       }
    848       if (havePS1) {
    849         secfilt[Nsec].flags |= ID_SECF_HAS_PS1;
    850       }
    851       if (haveUbercal) {
    852         secfilt[Nsec].flags |= ID_SECF_USE_UBERCAL;
    853       }
    854     }
    855   }
    856 
    857   if (isSetMrelFinal) {
    858     if (pass == 0) {
    859       DVOAverageFlags flagBits = ID_OBJ_EXT | ID_OBJ_EXT_ALT | ID_OBJ_GOOD | ID_OBJ_GOOD_ALT;
    860       // we attempt to set a few flags here; reset those bits before trying:
    861       average[0].flags &= ~flagBits;
    862     }
    863 
    864     if (NextPS1 && (NextPS1 > NpsfPS1)) {
    865       average[0].flags |= ID_OBJ_EXT;
    866     }
    867     if (GoodPS1) {
    868       average[0].flags |= ID_OBJ_GOOD;
    869     }
    870     if (Galaxy2MASS) {
    871       average[0].flags |= ID_OBJ_EXT_ALT;
    872     }
    873     if (Good2MASS) {
    874       average[0].flags |= ID_OBJ_GOOD_ALT;
    875     }
    876   }
    877   return (TRUE);
    878 }
    879 
    880753// only apply Stack operation on setMrelFinal in first pass
    881754// this function has 3 goals, not to be confused:
    … …  
    13381211  return (TRUE);
    13391212}
     1213
     1214int magStatsByRanking (StatDataSet *dataset, StatType *stats) {
     1215
     1216  liststats_init (stats);
     1217
     1218  if (dataset->Nlist == 0) return 0;
     1219
     1220  // we have a list of measurements for this Nsec value, along with errors, weights, and their rank
     1221  // first, sort by the rank (increasing)
     1222  sort_measures_by_ranking (dataset->ranking, dataset->measSeq, dataset->flxlist, dataset->errlist, dataset->wgtlist, dataset->Nlist);
     1223
     1224  // find the values with the same minimum rank:
     1225  int Nranking = 0;
     1226  int minRank = dataset->ranking[0]; // MIN (5, result->psfData[Nsec].ranking[0]);  -- only allow rank 5 or <
     1227
     1228  int i;
     1229  for (i = 0; (i < dataset->Nlist) && (dataset->ranking[i] == minRank); i++, Nranking++);
     1230
     1231  // if no valid measurements meet the minimum criteria, go to the next entry
     1232  if (Nranking < Nminmeas) {
     1233    return 0;
     1234  }
     1235
     1236  liststats (dataset->flxlist, dataset->errlist, dataset->wgtlist, Nranking, stats);
     1237  return (Nranking);
     1238}

  • branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelFinal.c

    r37037 r38293  
    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)
    … …  
    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  setMeasurementRank (catalog);
     105  setMrelOutput (catalog);
    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        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 (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 (measure[i].psfQFperf < 0.85)                 { measureRank[i] = 1; continue; }
     341       
     342    // RANK 0 : perfect measurement:
     343    measureRank[i] = 0;
     344  }
     345}
Note: See TracChangeset for help on using the changeset viewer.