IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38317


Ignore:
Timestamp:
May 25, 2015, 12:41:40 PM (11 years ago)
Author:
eugene
Message:

fixing various relphot issues

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

Legend:

Unmodified
Added
Removed

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

    r38316 r38317  
    106106    results->needSYN[i] = FALSE;
    107107    results->measSYN[i] = -1;
    108     results->minSYN[i]  = NAN;
    109108
    110109    results->psfQfMax[i]  = NAN;
    … …  
    151150    ALLOCATE (results->minSYN,  float, results->Nsecfilt);
    152151
    153     ALLOCATE (results->psfQfMax, float, results->Nsecfilt);
     152    ALLOCATE (results->psfQfMax,     float, results->Nsecfilt);
    154153    ALLOCATE (results->psfQfPerfMax, float, results->Nsecfilt);
    155154
    … …  
    158157    ALLOCATE (results->Next,           int, results->Nsecfilt);
    159158    ALLOCATE (results->NexpPS1,        int, results->Nsecfilt);
    160     ALLOCATE (results->havePS1,        int, results->Nsecfilt);
    161159    ALLOCATE (results->haveUbercal,    int, results->Nsecfilt);
    162160
    … …  
    199197    results->Next = NULL;
    200198    results->NexpPS1 = NULL;
    201     results->havePS1 = NULL;
    202199    results->haveUbercal = NULL;
    203200
    … …  
    231228  FREE (results->Next);
    232229  FREE (results->NexpPS1);
    233   FREE (results->havePS1);
    234230  FREE (results->haveUbercal);
    235231

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

    r38298 r38317  
    331331  }
    332332
    333   // if -synthphot_zpts is supplied, the map of corrections needed for the synthetic photometry is supplied
    334   // and used to tie down the synthetic magnitudes
    335   SYNTH_ZERO_POINTS = NULL;
    336   if ((N = get_argument (argc, argv, "-synthphot_zpts"))) {
    337     remove_argument (N, &argc, argv);
    338     myAssert (N < argc, "missing argument to -synthphot_means");
    339     SYNTH_ZERO_POINTS = strcreate (argv[N]);
    340     remove_argument (N, &argc, argv);
    341   }
    342 
    343333  // if -synthphot is chosen, an artificial w-band measurement is generated for each object
    344334  // with r and i photometry.  this was used to apply the ubercal photometry to w-band images.
    … …  
    384374    mode = UPDATE_AVERAGES;
    385375  }
     376
     377  // if -synthphot_zpts is supplied, the map of corrections needed for the synthetic photometry is supplied
     378  // and used to tie down the synthetic magnitudes
     379  // NOTE: SYNTH_ZERO_POINTS is used for the same file in both -synthphot_zpts and -synthphot_means
    386380  SYNTH_ZERO_POINTS = NULL;
     381  if ((N = get_argument (argc, argv, "-synthphot_zpts"))) {
     382    remove_argument (N, &argc, argv);
     383    myAssert (N < argc, "missing argument to -synthphot_zpts");
     384    SYNTH_ZERO_POINTS = strcreate (argv[N]);
     385    remove_argument (N, &argc, argv);
     386  }
    387387  if ((N = get_argument (argc, argv, "-synthphot_means"))) {
    388388    mode = SYNTH_PHOT;
    … …  
    495495  }
    496496  SYNTH_ZERO_POINTS = NULL;
     497  if ((N = get_argument (argc, argv, "-synthphot_zpts"))) {
     498    remove_argument (N, &argc, argv);
     499    myAssert (N < argc, "missing argument to -synthphot_zpts");
     500    SYNTH_ZERO_POINTS = strcreate (argv[N]);
     501    remove_argument (N, &argc, argv);
     502  }
    497503  if ((N = get_argument (argc, argv, "-synthphot_means"))) {
    498504    MODE = MODE_SYNTH_PHOT;

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

    r38314 r38317  
    6565      int Ns;
    6666
    67 # if (0)     
    68       DVOSecfiltFlags secfiltBits =
    69         ID_SECF_STAR_FEW |
    70         ID_SECF_STAR_POOR |
    71         ID_PHOTOM_PASS_0 |
    72         ID_PHOTOM_PASS_1 |
    73         ID_PHOTOM_PASS_2 |
    74         ID_PHOTOM_PASS_3 |
    75         ID_PHOTOM_PASS_4 |
    76         ID_SECF_USE_SYNTH |
    77         ID_SECF_USE_UBERCAL |
    78         ID_SECF_OBJ_EXT;
    79 # endif
    80 
    8167      for (Ns = 0; Ns < Nphotcodes; Ns++) {
    8268
    … …  
    8470        int Nsec = GetPhotcodeNsec(thisCode);
    8571
    86         // fprintf (stderr, "really use dvo_secfilt_init?");
    87         dvo_secfilt_init (&subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec]);
    88 
    89         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M      = NAN;
    90         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Map    = NAN;
    91         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].dM     = NAN;
    92         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Mstdev = NAN_S_SHORT;
    93         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Mchisq = NAN;
    94         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_20   = NAN_S_SHORT;
    95         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_80   = NAN_S_SHORT;
    96         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Ncode  = 0;
    97         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Nused  = 0;
    98         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~secfiltBits;
    99         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].ubercalDist = 1000;
     72        dvo_secfilt_init (&subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec], SECFILT_RESET_ALL);
    10073      }
    10174    }

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

    r35416 r38317  
    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 

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

    r38293 r38317  
    3737
    3838  myAssert (stats->statmode != STATS_NONE, "programming error, liststats mode not set");
     39
     40  liststats_init (stats);
    3941
    4042  ke = ks = dMo = 0;

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

    r38298 r38317  
    170170    free (BOUNDARY_TREE);
    171171    BOUNDARY_TREE = tmppath;
     172  }
     173  if (SYNTH_ZERO_POINTS) {
     174    char *tmppath = abspath(SYNTH_ZERO_POINTS, DVO_MAX_PATH);
     175    free (SYNTH_ZERO_POINTS);
     176    SYNTH_ZERO_POINTS = tmppath;
    172177  }
    173178
    … …  
    224229   
    225230    char tmpline[1024];
    226     if (VERBOSE)          { snprintf (tmpline, 1024, "%s -v",                   command);                                         strcpy (command, tmpline); }
    227     if (VERBOSE2)         { snprintf (tmpline, 1024, "%s -vv",                  command);                                         strcpy (command, tmpline); }
    228     if (RESET)            { snprintf (tmpline, 1024, "%s -reset",               command);                                         strcpy (command, tmpline); }
    229     if (RESET_ZEROPTS)    { snprintf (tmpline, 1024, "%s -reset-zpts",          command);                                         strcpy (command, tmpline); }
    230     if (UPDATE)           { snprintf (tmpline, 1024, "%s -update",              command);                                         strcpy (command, tmpline); }
    231     if (IS_DIFF_DB)       { snprintf (tmpline, 1024, "%s -is-diff-db",          command);                                         strcpy (command, tmpline); }
    232     if (!KEEP_UBERCAL)    { snprintf (tmpline, 1024, "%s -reset-ubercal",       command);                                         strcpy (command, tmpline); }
    233     if (UPDATE_CATFORMAT) { snprintf (tmpline, 1024, "%s -update-catformat %s", command, UPDATE_CATFORMAT);                       strcpy (command, tmpline); }
    234     if (BOUNDARY_TREE)    { snprintf (tmpline, 1024, "%s -boundary-tree %s",    command, BOUNDARY_TREE);                          strcpy (command, tmpline); }
     231    if (VERBOSE)           { snprintf (tmpline, 1024, "%s -v",                  command);                                         strcpy (command, tmpline); }
     232    if (VERBOSE2)          { snprintf (tmpline, 1024, "%s -vv",                 command);                                         strcpy (command, tmpline); }
     233    if (RESET)             { snprintf (tmpline, 1024, "%s -reset",              command);                                         strcpy (command, tmpline); }
     234    if (RESET_ZEROPTS)     { snprintf (tmpline, 1024, "%s -reset-zpts",         command);                                         strcpy (command, tmpline); }
     235    if (UPDATE)            { snprintf (tmpline, 1024, "%s -update",             command);                                         strcpy (command, tmpline); }
     236    if (IS_DIFF_DB)        { snprintf (tmpline, 1024, "%s -is-diff-db",         command);                                         strcpy (command, tmpline); }
     237    if (!KEEP_UBERCAL)     { snprintf (tmpline, 1024, "%s -reset-ubercal",      command);                                         strcpy (command, tmpline); }
     238    if (UPDATE_CATFORMAT)  { snprintf (tmpline, 1024, "%s -update-catformat %s", command, UPDATE_CATFORMAT);                      strcpy (command, tmpline); }
     239    if (BOUNDARY_TREE)     { snprintf (tmpline, 1024, "%s -boundary-tree %s",    command, BOUNDARY_TREE);                         strcpy (command, tmpline); }
    235240    if (SYNTH_ZERO_POINTS) { snprintf (tmpline, 1024, "%s -synthphot-zpts %s", command, SYNTH_ZERO_POINTS); strcpy (command, tmpline); }
    236241    if (SET_MREL_VERSION != 1) { snprintf (tmpline, 1024, "%s -set-mrel-version %d", command, SET_MREL_VERSION);                  strcpy (command, tmpline); }

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

    r38298 r38317  
    6262    if (RESET) {
    6363      Nsecfilt = catalog.Nsecfilt;
    64       DVOSecfiltFlags secfiltBits =
    65         ID_SECF_STAR_FEW    |
    66         ID_SECF_STAR_POOR   |
    67         ID_PHOTOM_PASS_0    |
    68         ID_PHOTOM_PASS_1    |
    69         ID_PHOTOM_PASS_2    |
    70         ID_PHOTOM_PASS_3    |
    71         ID_PHOTOM_PASS_4    |
    72         ID_SECF_USE_SYNTH   |
    73         ID_SECF_USE_UBERCAL |
    74         ID_SECF_HAS_STACK   |
    75         ID_SECF_HAS_PS1     |
    76         ID_SECF_OBJ_EXT     ;
    77      
     64
    7865      DVOAverageFlags astromBits =
    7966        ID_STAR_FIT_AVE   |
    … …  
    9178        catalog.average[j].stargal   = NAN;     // reset (will be re-calculated here)
    9279        for (k = 0; k < Nsecfilt; k++) {
    93           dvo_secfilt_init (&catalog.secfilt[j*Nsecfilt + k]);
    94           catalog.secfilt[j*Nsecfilt + k].flags &= ~secfiltBits; // XXX: only reset photom bits:
     80          dvo_secfilt_init (&catalog.secfilt[j*Nsecfilt + k], SECFILT_RESET_ALL);
    9581        }
    9682      }
    … …  
    10187    initMrel (&catalog, 1);
    10288    setMrelFinal (&catalog, NULL, TRUE);
     89    // XXX if we want to have options for setting warp, chip, stack independently, we need to init only the desired ones
    10390
    10491    if (!UPDATE) {

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

    r38316 r38317  
    1212    continue; }
    1313
    14 # define CHECK_VALID_MAG(MAG) (!isnan(MAG) && (MAG > 0.0) && (MAG < 30.0))
    15 # define CHECK_VALID_FLUX(MAG) (!isnan(FLUX))
    16 
    17 # if (0)
     14# define CHECK_VALID_MAG(MAG) (!isnan(MAG) && (MAG > -5.0) && (MAG < 30.0))
     15# define CHECK_VALID_FLUX(FLUX) (!isnan(FLUX))
     16
    1817static float MagToFlux (float Mag) {
    1918  float Flux = pow(10.0, -0.4*(Mag));
    2019  return (Flux);
    2120}
    22 # endif
    2321
    2422int print_measure_set_alt (Average *average, SecFilt *secfilt, Measure *measure) {
    … …  
    123121  int Galaxy2MASS = FALSE;
    124122  int haveTYCHO = FALSE;
     123
     124  float stargalmax = 0.0;
    125125
    126126  // assign measurements to the photcode lists
    … …  
    154154    results->Nmeas[Nsec] ++;
    155155
    156     // force the use of SYN even if we have PS1 mags?
    157     if (isSetMrelFinal && isGPC1synth(measureT[k].photcode)) {
    158       results->haveSYN[Nsec] = TRUE;
    159       results->measSYN[Nsec] = k;
    160       if (measureT[k].M < results->minSYN[Nsec]) results->needSYN[Nsec] = TRUE;
    161       continue;
    162     }
    163 
    164     // assign the max psfQF values for each Nsec value
     156    // various steps only relevant to the final pass
    165157    if (isSetMrelFinal) {
    166       if (measure[k].psfQF     > results->psfQfMax[Nsec])     results->psfQfMax[Nsec]     = measure[k].psfQF;
    167       if (measure[k].psfQFperf > results->psfQfPerfMax[Nsec]) results->psfQfPerfMax[Nsec] = measure[k].psfQFperf;
    168     }
    169 
    170     // are we a PS1 exposure photcode? (hard-wired photcodes)
    171     if (isGPC1chip(measureT[k].photcode)) { results->NexpPS1[Nsec] ++; results->havePS1[Nsec] = TRUE; }
     158      if (isfinite(measure[k].psfQF)     && (measure[k].psfQF     > results->psfQfMax[Nsec]))     results->psfQfMax[Nsec]     = measure[k].psfQF;
     159      if (isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > results->psfQfPerfMax[Nsec])) results->psfQfPerfMax[Nsec] = measure[k].psfQFperf;
     160
     161      // are we a PS1 exposure photcode? (hard-wired photcodes)
     162      if (isGPC1chip(measure[k].photcode)) {
     163        results->NexpPS1[Nsec] ++;
     164        results->havePS1[Nsec] = TRUE;
     165       
     166        // only count psfQF > 0.85 here
     167        if (isfinite(measure[k].extNsigma) && isfinite(measure[k].psfQF) && (measure[k].psfQF > 0.85)) {
     168          stargalmax = MAX (stargalmax, measure[k].extNsigma);
     169        }
     170      }
     171
     172      // force the use of SYN even if we have PS1 mags?
     173      if (isGPC1synth(measure[k].photcode)) {
     174        results->haveSYN[Nsec] = TRUE;
     175        results->measSYN[Nsec] = k;
     176        measureT[k].dbFlags |= ID_MEAS_SYNTH_MAG; // redundant with photcode
     177        // note that synthetic mags have the real mags in measure.M (no zero point offset)
     178        if (measure[k].M < results->minSYN[Nsec]) results->needSYN[Nsec] = TRUE;
     179        continue;
     180      }
     181    }
    172182
    173183    // ** Choose the calibration (depends on the mode : do I have an image reference or not?)
    … …  
    190200    }
    191201
    192     // XXX ubercal dist should be by Nsec
    193202    int myUbercalDist = getUbercalDist(meas, cat);
    194203    results->minUbercalDist[Nsec] = MIN(results->minUbercalDist[Nsec], myUbercalDist);
    … …  
    200209    int useUbercalWeight = isUbercal || (refPhotcode && (code->code == refPhotcode->code));
    201210
    202     float Map, dMap;
     211    float Map  = NAN;
     212    float dMap = NAN;
    203213    if (isSetMrelFinal) {
    204214      Map = PhotCat (&measure[k], MAG_CLASS_APER);
    … …  
    243253
    244254    // count the extended detections (all PS1 bands)
    245     if (isGPC1chip(measureT[k].photcode) && !isnan(Map) && !isnan(Mpsf)) {
     255    if (isSetMrelFinal && isGPC1chip(measureT[k].photcode) && !isnan(Map) && !isnan(Mpsf)) {
    246256      float dMagAp = Mpsf - Map;
    247257      float SigmaAp = hypot(0.1, 2.5*dMpsf);
    … …  
    260270  SynthZeroPoints *synthzpts = SynthZeroPointsGet();
    261271
     272  // find max values across filters
     273  float psfQfMax     = 0.0;
     274  float psfQfPerfMax = 0.0;
     275
    262276  // now calculate the mean stats for the Nsec bands.
    263277  for (Nsec = 0; Nsec < Nsecfilt; Nsec++) {
    264     dvo_secfilt_init (&secfilt[Nsec]); // XXX does this reset astrometry bits? no, right?
     278    dvo_secfilt_init (&secfilt[Nsec], SECFILT_RESET_CHIP); // this does not reset astrometry or STACK bits
    265279
    266280    // XXX hardwired
    … …  
    272286
    273287    if (haveTYCHO) {
    274       // secfilt[Nsec].flags |= ID_SECF_HAVE_TYCHO;
    275     }
    276 
    277     // XXX need to also supply ZP from map, apply tycho logic
     288      secfilt[Nsec].flags |= ID_SECF_HAS_TYCHO;
     289    }
     290
     291    // force the synthetic magnitude (apply zpt calibration if no tycho photometry exists)
    278292    if (results->needSYN[Nsec]) {
    279293      if (!isSetMrelFinal) continue;
    … …  
    300314        float *value = (float *) synthzpts->matrix[Nsec].buffer;
    301315        ZP = !isnan(value[Npix]) ? value[Npix] : 0.0;
     316        if (isfinite(value[Npix])) {
     317          secfilt[Nsec].flags |= ID_SECF_FIX_SYNTH;
     318        }
    302319      }
    303320      secfilt[Nsec].M      = Mpsf + ZP;
    … …  
    343360      secfilt[Nsec].psfQfPerfMax = results->psfQfPerfMax[Nsec];
    344361
     362      if (isfinite(secfilt[Nsec].psfQfMax)     && (secfilt[Nsec].psfQfMax     > psfQfMax))     psfQfMax     = secfilt[Nsec].psfQfMax;
     363      if (isfinite(secfilt[Nsec].psfQfPerfMax) && (secfilt[Nsec].psfQfPerfMax > psfQfPerfMax)) psfQfPerfMax = secfilt[Nsec].psfQfPerfMax;
     364
    345365      StatType *apstats = &results->apstats;
    346366      Nranking = magStatsByRanking (&results->aperData[Nsec], apstats);
    … …  
    370390      switch (minRank) {
    371391        case 0:
    372           secfilt[Nsec].flags |= ID_PHOTOM_PASS_0;
     392          secfilt[Nsec].flags |= ID_SECF_RANK_0;
    373393          if (results->havePS1[Nsec]) GoodPS1 = TRUE;
    374394          break;
    375395        case 1:
    376           secfilt[Nsec].flags |= ID_PHOTOM_PASS_1;
     396          secfilt[Nsec].flags |= ID_SECF_RANK_1;
    377397          if (results->havePS1[Nsec]) GoodPS1 = TRUE;
    378398          break;
    379399        case 2:
    380           secfilt[Nsec].flags |= ID_PHOTOM_PASS_2;
     400          secfilt[Nsec].flags |= ID_SECF_RANK_2;
    381401          if (results->havePS1[Nsec]) GoodPS1 = TRUE;
    382402          break;
    383403        case 3:
    384           secfilt[Nsec].flags |= ID_PHOTOM_PASS_3;
     404          secfilt[Nsec].flags |= ID_SECF_RANK_3;
    385405          break;
    386406        case 4:
    387           secfilt[Nsec].flags |= ID_PHOTOM_PASS_4;
     407          secfilt[Nsec].flags |= ID_SECF_RANK_4;
    388408          break;
    389409        default:
    … …  
    412432      average[0].flags |= ID_OBJ_GOOD_ALT;
    413433    }
     434    average[0].psfQF     = psfQfMax;
     435    average[0].psfQFperf = psfQfPerfMax;
    414436  }
    415437  return (TRUE);
    … …  
    636658  }
    637659
    638   setGlobalObjStats (average, measure);
    639 
    640660  return (TRUE);
    641661}
    … …  
    647667int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results) {
    648668
     669  // we are guaranteed to have average, measure, secfilt
     670  Average *average    = &catalog[0].average[ave];
    649671  off_t measureOffset = average->measureOffset;
    650672
    651   // we are guaranteed to have average, measure, secfilt
    652   Average *average  = &catalog[0].average[ave];
    653   Measure *measure  = &catalog[0].measure[measureOffset];
    654   SecFilt *secfilt  = &catalog[0].secfilt[ave*Nsecfilt];
    655   char *measureRank = &catalog[0].measureRank[measureOffset];
    656 
    657   off_t measureOffset = average->measureOffset;
    658 
    659   off_t k;
     673  Measure *measure    = &catalog[0].measure[measureOffset];
     674  SecFilt *secfilt    = &catalog[0].secfilt[ave*Nsecfilt];
     675  char *measureRank   = &catalog[0].measureRank[measureOffset];
     676
     677  off_t k, Nsec;
    660678  float Mcal = 0;
    661679
    … …  
    674692  for (Nsec = 0; Nsec < Nsecfilt; Nsec++) {
    675693    if (!IS_DIFF_DB) {
    676       myAssert (secfilt[Nsec].stackBestOff <= Nmeasure, "stackBestOff out of range");
    677       if (!FindImageSkycellID (secfilt[Nsec].stackBestOff, cat, &stkTessID[Nsec], &stkProjID[Nsec], &stkSkycellID[Nsec])) continue;
     694      myAssert (secfilt[Nsec].stackBestOff <= average->Nmeasure, "stackBestOff out of range");
     695      if (!FindImageSkycellID (secfilt[Nsec].stackBestOff, cat, &tessID, &projID, &skycellID)) continue;
    678696    }
    679697    results->tessID[Nsec]    = tessID;
    … …  
    684702  // assign measurements to the photcode lists
    685703  for (k = 0; k < average[0].Nmeasure; k++) {
     704    off_t meas = measureOffset + k;
    686705
    687706    // skip measurements that do not have a valid photcode (raise exception?)
    … …  
    702721
    703722    // use primary skycell for DIFF, stack skycell for WARP
    704     if (IS_DIFF_DB) {
    705       if (MatchImageSkycellID (meas, cat, priTessID, priProjID, priSkycellID)) {
    706         measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY;
    707       } else {
    708         continue;
    709       }
    710     } else {
    711       if (!MatchImageSkycellID (meas, cat, stkTessID, stkProjID, stkSkycellID)) continue;
     723    if (MatchImageSkycellID (meas, cat, results->tessID[Nsec], results->projID[Nsec], results->skycellID[Nsec])) {
     724      continue;
    712725    }
    713726
    714727    // assign the max psfQF values for each Nsec value
    715728    if (IS_DIFF_DB) {
     729      measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY;
    716730      if (measure[k].psfQF     > results->psfQfMax[Nsec])     results->psfQfMax[Nsec]     = measure[k].psfQF;
    717731      if (measure[k].psfQFperf > results->psfQfPerfMax[Nsec]) results->psfQfPerfMax[Nsec] = measure[k].psfQFperf;
    718732    }
    719733
    720     off_t meas = measureOffset + k;
    721734    if (getImageEntry (meas, cat) < 0) {
    722735      // measurements without an image are either external reference photometry or
    … …  
    731744    }
    732745    float Fcal = MagToFlux(-Mcal);
    733 
    734     // XX measure[k].dbFlags &= ~ID_MEAS_WARP_USED;
    735     // XXX measure[k].dbFlags |= ID_MEAS_WARP_USED;
    736746
    737747    // NOTE : I am using PhotCat not PhotSys for now since GPC1 chip-to-chip color terms
    … …  
    765775    if (CHECK_VALID_FLUX(Fkron)) {
    766776      int Nkron = results->kronData[Nsec].Nlist;
    767       results->kronData[Nsec].flxlist[Nkron] = Fap * Fcal;
    768       results->kronData[Nsec].errlist[Nkron] = dFap * Fcal;
     777      results->kronData[Nsec].flxlist[Nkron] = Fkron * Fcal;
     778      results->kronData[Nsec].errlist[Nkron] = dFkron * Fcal;
    769779      results->kronData[Nsec].wgtlist[Nkron] = 1.0;
    770780      results->kronData[Nsec].ranking[Nkron] = measureRank ? measureRank[k] : 0;
    … …  
    774784  }
    775785
     786  // find max values across filters
     787  float psfQfMax     = 0.0;
     788  float psfQfPerfMax = 0.0;
     789
    776790  // now calculate the mean stats for the Nsec bands.
    777791  for (Nsec = 0; Nsec < Nsecfilt; Nsec++) {
    778     // XXX cannot do this -- resets psf mags
    779     // dvo_secfilt_init (&secfilt[Nsec]); // XXX does this reset astrometry bits? no, right?
     792    dvo_secfilt_init (&secfilt[Nsec], SECFILT_RESET_WARP);
    780793
    781794    if (IS_DIFF_DB) {
    … …  
    783796      secfilt[Nsec].psfQfMax     = results->psfQfMax[Nsec];
    784797      secfilt[Nsec].psfQfPerfMax = results->psfQfPerfMax[Nsec];
     798
     799      if (isfinite(secfilt[Nsec].psfQfMax)     && (secfilt[Nsec].psfQfMax     > psfQfMax))     psfQfMax     = secfilt[Nsec].psfQfMax;
     800      if (isfinite(secfilt[Nsec].psfQfPerfMax) && (secfilt[Nsec].psfQfPerfMax > psfQfPerfMax)) psfQfPerfMax = secfilt[Nsec].psfQfPerfMax;
    785801    }
    786802
    … …  
    799815
    800816    // if too few valid measurements meet the minimum criteria, go to the next entry
    801     StatType *aperstats = &results->aperstats;
    802     Nranking = magStatsByRanking (&results->aperData[Nsec], aperstats);
     817    StatType *apstats = &results->apstats;
     818    Nranking = magStatsByRanking (&results->aperData[Nsec], apstats);
    803819    if (Nranking) {
    804       secfilt[Nsec].FapWrp     = aperstats->mean;
    805       secfilt[Nsec].dFapWrp    = aperstats->error;
    806       secfilt[Nsec].sFapWrp    = (aperstats->Nmeas > 1) ? aperstats->chisq : NAN;
    807       secfilt[Nsec].NusedApWrp = aperstats->Nmeas;
     820      secfilt[Nsec].FapWrp     = apstats->mean;
     821      secfilt[Nsec].dFapWrp    = apstats->error;
     822      secfilt[Nsec].sFapWrp    = (apstats->Nmeas > 1) ? apstats->chisq : NAN;
     823      secfilt[Nsec].NusedApWrp = apstats->Nmeas;
    808824      secfilt[Nsec].MapWrp     = isnan(secfilt[Nsec].FapWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FapWrp); // 8.9 since flux is in Jy
    809825    }
    … …  
    823839    secfilt[Nsec].NwarpGood = results->NmeasGood[Nsec];
    824840  }
    825 
    826   setGlobalObjStats (average, measure);
    827 
    828   return (TRUE);
    829 }
    830 
    831 int setGlobalObjStats (Average *average, Measure *measure) {
    832 
    833   off_t k;
    834 
    835   myAbort ("this function makes no sense");
    836 
    837   float psfQFmax     = 0.0;
    838   float psfQFperfmax = 0.0;
    839   float stargalmax   = 0.0;
    840 
    841   for (k = 0; k < average[0].Nmeasure; k++) {
    842 
    843     if (isfinite(measure[k].psfQF)) {
    844       psfQFmax     = MAX (psfQFmax,     measure[k].psfQF);
    845     }
    846     if (isfinite(measure[k].psfQFperf)) {
    847       psfQFperfmax = MAX (psfQFperfmax, measure[k].psfQFperf);
    848     }
    849 
    850     // only count psfQF > 0.85 here
    851     if (isfinite(measure[k].extNsigma) && (measure[k].psfQF > 0.85)) {
    852       stargalmax   = MAX (stargalmax,   measure[k].extNsigma);
    853     }
    854   }
    855   average[0].psfQF     = psfQFmax;
    856   average[0].psfQFperf = psfQFperfmax;
    857   average[0].stargal   = stargalmax;
    858 
    859841  return (TRUE);
    860842}

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

    r38297 r38317  
    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;
Note: See TracChangeset for help on using the changeset viewer.