IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38314


Ignore:
Timestamp:
May 24, 2015, 10:18:27 AM (11 years ago)
Author:
eugene
Message:

mark used measurements; re-work the warp analysis

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

Legend:

Unmodified
Added
Removed

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

    r38298 r38314  
    137137
    138138  int *Nmeas;                 // count of PS1 exposure (chip) measurements for this secfilt
     139  int *NmeasGood;             // count of PS1 exposure (chip) measurements for this secfilt
    139140  int *Next;                  // count of PS1 exposure (chip) measurements for this secfilt
    140141  int *NexpPS1;               // count of PS1 exposure (chip) measurements for this secfilt

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

    r38293 r38314  
    692692      list[N] = Msys - Mrel - Mmos - Mgrid + Mflat;
    693693      dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
     694      if (fabs(list[N]) > 0.03) {
     695        // fprintf (stderr, "deviant\n");
     696      }
     697
    694698      if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
    695699        Mlist[Nbright] = list[N];
    … …  
    722726
    723727    if (PLOTSTUFF) {
    724       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);
    725729      plot_setMcal (list, N, &stats, CLOUD_TOLERANCE);
    726730    }

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

    r38298 r38314  
    112112
    113113    results->Nmeas[i] = 0;
     114    results->NmeasGood[i] = 0;
    114115    results->Next[i] = 0;
    115116    results->NexpPS1[i] = 0;
    … …  
    149150
    150151    ALLOCATE (results->Nmeas,          int, results->Nsecfilt);
     152    ALLOCATE (results->NmeasGood,      int, results->Nsecfilt);
    151153    ALLOCATE (results->Next,           int, results->Nsecfilt);
    152154    ALLOCATE (results->NexpPS1,        int, results->Nsecfilt);
    … …  
    12401242            }
    12411243            Ndel ++;
     1244          } else {
     1245            catalog[i].measureT[ilist[k]].dbFlags &= ~ID_MEAS_POOR_PHOTOM;
     1246            if (final) {
     1247              // for the final pass, we have a duplicate set of values in measure and measureT
     1248              catalog[i].measure[ilist[k]].dbFlags &= ~ID_MEAS_POOR_PHOTOM;
     1249            }
    12421250          }
    12431251        }

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

    r37116 r38314  
    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++) {

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

    r35759 r38314  
    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

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

    r38298 r38314  
    11# include "relphot.h"
     2int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags flags);
     3
    24# define UBERCAL_WEIGHT 100.0
    35void sort_StatDataSet (StatDataSet *dataset);
    … …  
    1113
    1214# define CHECK_VALID_MAG(MAG) (!isnan(MAG) && (MAG > 0.0) && (MAG < 30.0))
     15# define CHECK_VALID_FLUX(MAG) (!isnan(FLUX))
    1316
    1417# if (0)
    … …  
    318321
    319322    if (isSetMrelFinal) {
     323      if ((average->objID == 0x1) && (average->catID)) {
     324        fprintf (stderr, "test obj\n");
     325      }
     326
     327      // mark the measurements matching this ranking
     328      markMeasureByRanking (&results->psfData[Nsec], measure, minRank, ID_MEAS_PHOTOM_PSF);
     329
    320330      // when running -averages, we have no information about the images, so we cannot set this
    321331      if (results->minUbercalDist[Nsec] > -1) {
    … …  
    341351        secfilt[Nsec].NusedAp = Nranking;
    342352      }
     353      markMeasureByRanking (&results->aperData[Nsec], measure, minRank, ID_MEAS_PHOTOM_APER);
    343354
    344355      StatType *kronstats = &results->kronstats;
    … …  
    350361        secfilt[Nsec].NusedKron = Nranking;
    351362      }
     363      markMeasureByRanking (&results->kronData[Nsec], measure, minRank, ID_MEAS_PHOTOM_KRON);
    352364
    353365      // does this object appear extended in > 50% of measurements?
    … …  
    629641}
    630642
    631 # undef SKIP_THIS_MEAS
    632 # define SKIP_THIS_MEAS(REASON) {                               \
    633     measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM;  \
    634     results->REASON ++;                                         \
    635     continue; }
    636 
    637643// set mean of forced-warp measurements (selected by photcode range for now):
    638644// somewhat simplified relative to chip-photometry:
    … …  
    640646int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results) {
    641647
    642 # if (0)
    643 
    644648  off_t measureOffset = average->measureOffset;
    645649
    … …  
    652656  float Mcal = 0;
    653657
    654   // we are measuring means for 3 types of FLUXes: psf, ap, kron.  I am using the psf mag
    655   // error for the ap mags, but krons have their own errors.  it is an open question if I
    656   // should be doing weighted or unweighted fits (this is a user option)
    657   double *Fpsflist  = results->Mpsflist;
    658   double *dpsflist  = results->dpsflist;
    659   double *wpsflist  = results->wpsflist;
    660 
    661   double *Faplist   = results->Maplist;
    662   double *daplist   = results->daplist;
    663   double *waplist   = results->waplist;
    664 
    665   double *Fkronlist = results->Mkronlist;
    666   double *dkronlist = results->dkronlist;
    667   double *wkronlist = results->wkronlist;
    668 
    669   StatType *psfstats  = &results->psfstats;
    670   StatType *apstats   = &results->apstats;
    671   StatType *kronstats = &results->kronstats;
    672 
    673658  // option for a test print
    674659  if ((average[0].objID == 0xd6e) && (average[0].catID == 0x4984)) {
    … …  
    683668  }
    684669
    685   int Ns;
    686   for (Ns = 0; Ns < Nphotcodes; Ns++) {
    687 
    688     int thisCode = photcodes[Ns][0].code;
    689     int Nsec = GetPhotcodeNsec(thisCode);
    690 
    691     /* calculate the average mag in this SEC photcode for a single star */
    692 
    693     /* star/photcodes already calibrated */
    694     if (found[Nsec]) continue; 
    695      
     670  for (Nsec = 0; Nsec < Nsecfilt; Nsec++) {
    696671    // I need to restrict the measurements coming only from the same skycell as the secfilt values
    697672    // if we cannot make that association, skip this set of warps
    698673    int stkTessID, stkProjID, stkSkycellID;
    699 
     674   
    700675    if (!IS_DIFF_DB) {
    701676      myAssert (secfilt[Nsec].stackBestOff <= Nmeasure, "stackBestOff out of range");
    702       if (!FindImageSkycellID (secfilt[Nsec].stackBestOff, cat, &stkTessID, &stkProjID, &stkSkycellID)) continue;
    703     }
    704 
    705     float psfQfMax = 0.0;
    706     float psfQfPerfMax = 0.0;
    707 
    708     off_t meas = measureOffset;
    709 
    710     int Nap = 0;
    711     int Npsf = 0;
    712     int Nkron = 0;
    713 
    714     int Nwarp = 0;
    715     int NwarpGood = 0;
    716     for (k = 0; k < average[0].Nmeasure; k++, meas++) {
    717 
    718       // only examine gpc1 forced-warp data
    719       if (!isGPC1warp(measure[k].photcode)) continue;
    720 
    721       // skip measurements that do not match the current photcode
    722       PhotCode *code = GetPhotcodebyCode (measure[k].photcode);
    723       if (!code) continue;
    724       if (code->equiv != thisCode) { continue; }
    725 
    726       Nwarp ++;
    727       if (measure[k].psfQFperf > 0.85) NwarpGood ++;
    728 
    729       measure[k].dbFlags &= ~ID_MEAS_WARP_USED;
    730 
    731       if (IS_DIFF_DB) {
    732         if (MatchImageSkycellID (meas, cat, priTessID, priProjID, priSkycellID)) {
    733           measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY;
    734         } else {
    735           continue;
    736         }
     677      if (!FindImageSkycellID (secfilt[Nsec].stackBestOff, cat, &stkTessID[Nsec], &stkProjID[Nsec], &stkSkycellID[Nsec])) continue;
     678    }
     679  }
     680
     681  // assign measurements to the photcode lists
     682  for (k = 0; k < average[0].Nmeasure; k++) {
     683
     684    // skip measurements that do not have a valid photcode (raise exception?)
     685    PhotCode *code = GetPhotcodebyCode (measure[k].photcode);
     686    myAssert (code, "invalid photcode??");
     687
     688    // only examine gpc1 forced-warp data
     689    if (!isGPC1warp(measure[k].photcode)) continue;
     690
     691    // assign the Nsec value so we can assign to the right lists
     692    int Nsec = GetPhotcodeNsec (code->equiv);
     693    if (Nsec < 0) continue; // skip measurements which do not have an equiv average photcode
     694
     695    results->Nmeas[Nsec] ++;
     696    if (measure[k].psfQFperf > 0.85) results->NwarpGood[Nsec] ++;
     697
     698    measure[k].dbFlags &= ~ID_MEAS_WARP_USED;
     699
     700    // use primary skycell for DIFF, stack skycell for WARP
     701    if (IS_DIFF_DB) {
     702      if (MatchImageSkycellID (meas, cat, priTessID, priProjID, priSkycellID)) {
     703        measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY;
    737704      } else {
    738         if (!MatchImageSkycellID (meas, cat, stkTessID, stkProjID, stkSkycellID)) continue;
    739       }
    740  
    741       if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad);
    742 
    743       if (getImageEntry (meas, cat) < 0) {
    744         // measurements without an image are either external reference photometry or
    745         // data for which the associated image has not been loaded (probably because of
    746         // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
    747         Mcal = measure[k].Mcal; // check that this is zero for loaded REF value
    748       } else {
    749         // use getMcal not getMcal_alt?
    750         Mcal  = getMcal_alt (meas, cat, NULL, measure[k].Xccd, measure[k].Yccd);
    751         // Mcal  = getMcal (meas, cat);
    752         if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
    753       }
    754       float Fcal = MagToFlux(-Mcal);
    755 
    756       // these are only used below if IS_DIFF_DB
    757       if (IS_DIFF_DB) {
    758         if (measure[k].psfQF     > psfQfMax)     psfQfMax     = measure[k].psfQF;
    759         if (measure[k].psfQFperf > psfQfPerfMax) psfQfPerfMax = measure[k].psfQFperf;
    760       }
    761 
    762       // in the calculations below,
    763       // ...list gives the error per measurement, wlist gives the weight
    764       // we can modify the error and weight in a few ways:
    765       // 1) MIN_ERROR guarantees a floor
    766       // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
    767 
    768       // skip some absurd values NAN
    769       // NOTE : I am using PhotCat not PhotSys for now since GPC1 chip-to-chip color terms
    770       // are small (and not measured)
    771       float Fpsf = PhotFluxCat (&measure[k], MAG_CLASS_PSF);
    772       if (isnan(Fpsf)) SKIP_THIS_MEAS(Nsys);
    773       // if (Msys <  0.0) SKIP_THIS_MEAS(Nsys);
    774       // if (Msys > 30.0) SKIP_THIS_MEAS(Nsys);
    775 
    776       float dFpsf = PhotFluxCatErr (&measure[k], MAG_CLASS_PSF);
    777       if (isnan(dFpsf)) SKIP_THIS_MEAS (Nsys);
    778 
    779       int isBad = (measure[k].psfQF < 0.85);
    780       isBad |= isnan(measure[k].psfQF);
    781       if (isBad) SKIP_THIS_MEAS (Nsys);
    782 
    783       int isPoor = (measure[k].psfQFperf < 0.85);
    784       if ((pass == 0) && isPoor) SKIP_THIS_MEAS (Nsys);
    785 
    786       dFpsf = MAX (dFpsf, MIN_ERROR*Fpsf); // MIN_ERROR is a fractional error
    787       Fpsflist[Npsf] = Fpsf * Fcal;
    788       dpsflist[Npsf] = dFpsf * Fcal;
    789       wpsflist[Npsf] = 1.0;
    790       Npsf ++;
    791 
    792       measure[k].dbFlags |= ID_MEAS_WARP_USED;
    793 
    794       float Fap = PhotFluxCat (&measure[k], MAG_CLASS_APER);
    795       float dFap = PhotFluxCatErr (&measure[k], MAG_CLASS_APER);
    796       if (!isnan(Fap)) {
    797         Faplist[Nap] = Fap * Fcal;
    798         daplist[Nap] = dFap * Fcal;
    799         waplist[Nap] = 1.0; // drop weight lists?
    800         Nap ++;
    801       }
    802 
    803       float Fkron = PhotFluxCat (&measure[k], MAG_CLASS_KRON);
    804       float dFkron = PhotFluxCatErr (&measure[k], MAG_CLASS_KRON);
    805       if (!isnan(Fkron)) {
    806         Fkronlist[Nkron] = Fkron * Fcal;
    807         dkronlist[Nkron] = dFkron * Fcal;
    808         wkronlist[Nkron] = 1.0;
    809         Nkron ++;
    810       }
    811     }
    812     if (Npsf < 1) {
    813       secfilt[Nsec].Nwarp        = Nwarp;
    814       secfilt[Nsec].NwarpGood    = NwarpGood;
    815       secfilt[Nsec].NusedWrp     = Npsf;
    816       secfilt[Nsec].NusedApWrp   = Nap;
    817       secfilt[Nsec].NusedKronWrp = Nkron;
    818       continue;
    819     }
    820 
    821     found[Nsec] = TRUE;
    822 
    823     liststats (Fpsflist, dpsflist, wpsflist, Npsf, psfstats);
     705        continue;
     706      }
     707    } else {
     708      if (!MatchImageSkycellID (meas, cat, stkTessID, stkProjID, stkSkycellID)) continue;
     709    }
     710
     711    // assign the max psfQF values for each Nsec value
     712    if (IS_DIFF_DB) {
     713      if (measure[k].psfQF     > results->psfQfMax[Nsec])     results->psfQfMax[Nsec]     = measure[k].psfQF;
     714      if (measure[k].psfQFperf > results->psfQfPerfMax[Nsec]) results->psfQfPerfMax[Nsec] = measure[k].psfQFperf;
     715    }
     716
     717    if (getImageEntry (meas, cat) < 0) {
     718      // measurements without an image are either external reference photometry or
     719      // data for which the associated image has not been loaded (probably because of
     720      // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
     721      Mcal = measure[k].Mcal; // check that this is zero for loaded REF value
     722    } else {
     723      // use getMcal not getMcal_alt?
     724      Mcal  = getMcal_alt (meas, cat, NULL, measure[k].Xccd, measure[k].Yccd);
     725      // Mcal  = getMcal (meas, cat);
     726      if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
     727    }
     728    float Fcal = MagToFlux(-Mcal);
     729
     730    // XX measure[k].dbFlags &= ~ID_MEAS_WARP_USED;
     731    // XXX measure[k].dbFlags |= ID_MEAS_WARP_USED;
     732
     733    // NOTE : I am using PhotCat not PhotSys for now since GPC1 chip-to-chip color terms
     734    // are small (and not measured)
     735    float Fpsf = PhotFluxCat (&measure[k], MAG_CLASS_PSF);
     736    float dFpsf = PhotFluxCatErr (&measure[k], MAG_CLASS_PSF);
     737    if (CHECK_VALID_FLUX(Fpsf)) {
     738      int Npsf = results->psfData[Nsec].Nlist;
     739      results->psfData[Nsec].flxlist[Npsf] = Fpsf * Fcal;
     740      results->psfData[Nsec].errlist[Npsf] = dFpsf * Fcal;
     741      results->psfData[Nsec].wgtlist[Npsf] = 1.0;
     742      results->psfData[Nsec].ranking[Npsf] = measureRank ? measureRank[k] : 0;
     743      results->psfData[Nsec].measSeq[Npsf] = k;
     744      results->psfData[Nsec].Nlist ++;
     745    }     
     746
     747    float Fap = PhotFluxCat (&measure[k], MAG_CLASS_APER);
     748    float dFap = PhotFluxCatErr (&measure[k], MAG_CLASS_APER);
     749    if (CHECK_VALID_FLUX(Fap)) {
     750      int Naper = results->aperData[Nsec].Nlist;
     751      results->aperData[Nsec].flxlist[Naper] = Fap * Fcal;
     752      results->aperData[Nsec].errlist[Naper] = dFap * Fcal;
     753      results->aperData[Nsec].wgtlist[Naper] = 1.0;
     754      results->aperData[Nsec].ranking[Naper] = measureRank ? measureRank[k] : 0;
     755      results->aperData[Nsec].measSeq[Naper] = k;
     756      results->aperData[Nsec].Nlist ++;
     757    }
     758
     759    float Fkron = PhotFluxCat (&measure[k], MAG_CLASS_KRON);
     760    float dFkron = PhotFluxCatErr (&measure[k], MAG_CLASS_KRON);
     761    if (CHECK_VALID_FLUX(Fkron)) {
     762      int Nkron = results->kronData[Nsec].Nlist;
     763      results->kronData[Nsec].flxlist[Nkron] = Fap * Fcal;
     764      results->kronData[Nsec].errlist[Nkron] = dFap * Fcal;
     765      results->kronData[Nsec].wgtlist[Nkron] = 1.0;
     766      results->kronData[Nsec].ranking[Nkron] = measureRank ? measureRank[k] : 0;
     767      results->kronData[Nsec].measSeq[Nkron] = k;
     768      results->kronData[Nsec].Nlist ++;
     769    }
     770  }
     771
     772  // now calculate the mean stats for the Nsec bands.
     773  for (Nsec = 0; Nsec < Nsecfilt; Nsec++) {
     774    // XXX cannot do this -- resets psf mags
     775    // dvo_secfilt_init (&secfilt[Nsec]); // XXX does this reset astrometry bits? no, right?
    824776
    825777    if (IS_DIFF_DB) {
    826778      // for non DIFF_DB, these are set in setMrelAverageExposure
    827       secfilt[Nsec].psfQfMax     = psfQfMax;
    828       secfilt[Nsec].psfQfPerfMax = psfQfPerfMax;
    829     }
    830 
    831     secfilt[Nsec].FpsfWrp  = psfstats->mean;
    832     secfilt[Nsec].dFpsfWrp = psfstats->error;
    833     secfilt[Nsec].sFpsfWrp = psfstats->sigma; // Mstdev is in millimags (not enough space for more precision)
    834     secfilt[Nsec].NusedWrp = psfstats->Nmeas;
    835     secfilt[Nsec].MpsfWrp  = isnan(secfilt[Nsec].FpsfWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FpsfWrp); // 8.9 since flux is in Jy
    836 
    837     // NOTE : use the modified weight for apmags as well as psf mags
    838     liststats (Faplist, daplist, waplist, Nap, apstats);
    839     secfilt[Nsec].FapWrp     = Nap > 0 ? apstats->mean : NAN;
    840     secfilt[Nsec].dFapWrp    = Nap > 0 ? apstats->error : NAN;
    841     secfilt[Nsec].sFapWrp    = Nap > 0 ? apstats->sigma : NAN;
    842     secfilt[Nsec].NusedApWrp = Nap;
    843     secfilt[Nsec].MapWrp  = isnan(secfilt[Nsec].FapWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FapWrp); // 8.9 since flux is in Jy
    844 
    845     liststats (Fkronlist, dkronlist, wkronlist, Nkron, kronstats);
    846     secfilt[Nsec].FkronWrp     = Nkron > 0 ? kronstats->mean  : NAN;
    847     secfilt[Nsec].dFkronWrp    = Nkron > 0 ? kronstats->error : NAN;
    848     secfilt[Nsec].sFkronWrp    = Nkron > 0 ? kronstats->sigma : NAN;
    849     secfilt[Nsec].NusedKronWrp = Nkron;
    850     secfilt[Nsec].MkronWrp  = isnan(secfilt[Nsec].FkronWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FkronWrp); // 8.9 since flux is in Jy
    851 
    852     secfilt[Nsec].Nwarp     = Nwarp;
    853     secfilt[Nsec].NwarpGood = NwarpGood;
     779      secfilt[Nsec].psfQfMax     = results->psfQfMax[Nsec];
     780      secfilt[Nsec].psfQfPerfMax = results->psfQfPerfMax[Nsec];
     781    }
     782
     783    // if too few valid measurements meet the minimum criteria, go to the next entry
     784    StatType *psfstats = &results->psfstats;
     785    int Nranking = magStatsByRanking (&results->psfData[Nsec], psfstats);
     786    if (Nranking) {
     787      secfilt[Nsec].FpsfWrp  = psfstats->mean;
     788      secfilt[Nsec].dFpsfWrp = psfstats->error;
     789      secfilt[Nsec].sFpsfWrp = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
     790      secfilt[Nsec].NusedWrp = psfstats->Nmeas;
     791      secfilt[Nsec].MpsfWrp  = isnan(secfilt[Nsec].FpsfWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FpsfWrp); // 8.9 since flux is in Jy
     792    }
     793    int minRank = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10;
     794    markMeasureByRanking (&results->psfData[Nsec], measure, minRank, IID_MEAS_WARP_USED);
     795
     796    // if too few valid measurements meet the minimum criteria, go to the next entry
     797    StatType *aperstats = &results->aperstats;
     798    Nranking = magStatsByRanking (&results->aperData[Nsec], aperstats);
     799    if (Nranking) {
     800      secfilt[Nsec].FapWrp     = aperstats->mean;
     801      secfilt[Nsec].dFapWrp    = aperstats->error;
     802      secfilt[Nsec].sFapWrp    = (aperstats->Nmeas > 1) ? aperstats->chisq : NAN;
     803      secfilt[Nsec].NusedApWrp = aperstats->Nmeas;
     804      secfilt[Nsec].MapWrp     = isnan(secfilt[Nsec].FapWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FapWrp); // 8.9 since flux is in Jy
     805    }
     806
     807    // if too few valid measurements meet the minimum criteria, go to the next entry
     808    StatType *kronstats = &results->kronstats;
     809    Nranking = magStatsByRanking (&results->kronData[Nsec], kronstats);
     810    if (Nranking) {
     811      secfilt[Nsec].FkronWrp     = kronstats->mean;
     812      secfilt[Nsec].dFkronWrp    = kronstats->error;
     813      secfilt[Nsec].sFkronWrp    = (kronstats->Nmeas > 1) ? kronstats->chisq : NAN;
     814      secfilt[Nsec].NusedKronWrp = kronstats->Nmeas;
     815      secfilt[Nsec].MkronWrp     = isnan(secfilt[Nsec].FkronWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FkronWrp); // 8.9 since flux is in Jy
     816    }
     817
     818    secfilt[Nsec].Nwarp     = results->Nmeas[Nsec];
     819    secfilt[Nsec].NwarpGood = results->NmeasGood[Nsec];
    854820  }
    855821
    856822  setGlobalObjStats (average, measure);
    857 
    858 # endif
    859823
    860824  return (TRUE);
    … …  
    909873  liststats (dataset->flxlist, dataset->errlist, dataset->wgtlist, Nranking, stats);
    910874  return (Nranking);
     875}
     876
     877int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags flags) {
     878  int i;
     879
     880  for (i = 0; i < dataset->Nlist; i++) {
     881    int k = dataset->measSeq[i];
     882    if (dataset->ranking[i] > minrank) {
     883      measure[k].dbFlags &= ~flags;
     884    } else {
     885      measure[k].dbFlags |= flags;
     886    }
     887  }
     888  return TRUE;
    911889}
    912890
Note: See TracChangeset for help on using the changeset viewer.