IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 24982


Ignore:
Timestamp:
Aug 3, 2009, 4:08:06 AM (17 years ago)
Author:
Sebastian Jester
Message:

Only load objects marked OK_RUN in the SDSS tsObj table's status bitmask to avoid counting the same photons more than once. Copious comments (marked by XXXJester) with answers to previous questions in previous comments, new questions, and suggestions about various issues.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/sj_branches/sj_SDSSaddstar_branch_r23928_20090420/Ohana/src/addstar/src/ReadStarsSDSS.c

    r23943 r24982  
    6767  strcpy (filtname[4], "z");
    6868
    69   // XXXYYY Why these? and not the 'ref' version? see 'Why these?' below, too.
     69  // XXXJester: SDSS photometry is read with 'dep' photcodes - recalibratable in principle, not standard-star
     70  // photometry (which would get 'ref').
    7071  NAMED_PHOTCODE_AND_ZP (photcode[0], zeropt[0], "U_SDSS");
    7172  NAMED_PHOTCODE_AND_ZP (photcode[1], zeropt[1], "G_SDSS");
    … …  
    7475  NAMED_PHOTCODE_AND_ZP (photcode[4], zeropt[4], "Z_SDSS");
    7576
    76   // XXXYYYZZZ SDSS tables have special flags for undefined/unmeasured values and errors:
    77   // -9999 flags unmeasured values, and the corresponding error may or may not be meaningful
    78   // -1000 flags errors that are not determined, even though the corresponding quantity is.
    79   // These special values need to be trapped here to avoid averaging meaningful numbers with the flag values.
     77  // XXXJester:ZZZ SDSS tables have special flags for undefined/unmeasured values and errors: -9999 flags
     78  // unmeasured values, and the corresponding error may or may not be meaningful -1000 flags errors that are
     79  // not determined, even though the corresponding quantity is.  These special values need to be trapped here
     80  // to avoid averaging meaningful numbers with the flag values.  However, we can't just skip an entire star
     81  // just because one measurement is missing in one quantity. Perhaps the right thing to do is to replace
     82  // -9999 and -1000 by NaN; however, does the averaging routine (i.e. liststats() in
     83  // Ohana/src/relphot/src/liststats.c, or maybe rather the sort functions in there) correctly ignore NaNs? I don't
     84  // think so. In fact, I think not ignoring NaNs is a generic problem, because NaNs tend to be contagious,
     85  // and just a single bad/missing measurement should not invalidate the average quantity of anything.
    8086
    8187  // various header values needed to calculate per-star data below
    … …  
    94100  tzero[4] = ohana_mjd_to_sec (mjd[4]);
    95101
    96   // XXXYYY EDR paper says "psfWidth is is the effective width also determined at the center of each frame. It
     102  // XXXJester: EDR paper says "psfWidth is is the effective width also determined at the center of each frame. It
    97103  // is a good generic number to quote for the seeing on each frame." It is stored in the tsField file as
    98   // psf_width.
     104  // psf_width. If we read the actual per-field zero points and extinction coefficients, we could also read psf_width.
    99105  gfits_scan (table.header, "SEEING_U", "%f", 1, &seeing[0]);
    100106  gfits_scan (table.header, "SEEING_G", "%f", 1, &seeing[1]);
    … …  
    103109  gfits_scan (table.header, "SEEING_Z", "%f", 1, &seeing[4]);
    104110
    105   // XXXYYY What's the meaning of photErr? Where is it used again?
     111  // XXXJester: What's the meaning of photErr? Where is it used again?
    106112  gfits_scan (table.header, "PSFERR_U", "%f", 1, &photErr[0]);
    107113  gfits_scan (table.header, "PSFERR_G", "%f", 1, &photErr[1]);
    … …  
    112118  gfits_scan (header, "CAMCOL", "%d", 1, &camcol); // value in header is usec / unbinned row
    113119
    114   // XXXYYY which number is read here? In other words, where is the zeropoint set that is returned by
    115   // GetZeroPoint()? I can't find the corresponding SetZeroPoint()
     120  // XXXJester: The corresponding SetZeroPoint() is in ConfigInit.c, where ZERO_PT from ~/.ptolemyrc is read
    116121  ZeroPt = GetZeroPoint();
    117122
    … …  
    147152  // compensate for the difference.
    148153
    149   Nstars_used = Nstars;
     154  // XXXYYYJester: the actual C_0 and the airmass coefficient are reported in the tsField file that go with
     155  // every tsObj file (i.e. C_0 = - aa - kk * airmass). Worse, SDSS reports *luptitudes* (asinh magnitudes)
     156  // not conventional ("Pogson") magnitudes, so in fact SDSS reports -(2.5/ln(10))*[asinh((f/f0)/2b)+ln(b)]
     157  // where f/f0 = counts/sec * 10^[0.4*(aa + kk * airmass)]. See
     158  // http://www.sdss.org/dr7/algorithms/fluxcal.html#counts2mag
     159
     160  Nstars_used = 0;
    150161  for (i = 0; i < Nstars; i++) {
    151     // Skip objects that are not OK_RUN
    152     if ( ((status[i]) & AR_OBJECT_STATUS_OK_RUN) != 0) {
    153       Nstars_used--;
     162    // Skip objects that are not OK_RUN (including other objects leads to counting the same photons more than
     163    // once - see comment near images[N].NY = ... below)
     164    if ( ((status[i]) & AR_OBJECT_STATUS_OK_RUN) == 0) {
    154165      continue;
    155166    }
     167    Nstars_used += 1;
    156168    for (j = 0; j < NFILTER; j++) {
    157       N = NFILTER*i + j;
     169      N = NFILTER * (Nstars_used-1) + j;
    158170      InitStar (&stars[N]);
    159171
    160172      // any values not explicitly set are left at 0.0
    161       stars[N].average.R         = ra[i] + dCOS(dec[i]) * offsetRa[N] / 3600.0;
     173      // XXXYYYJester: the offsetRa here used to be multiplied by dCOS(dec[i]), but that's already been done
     174      // by SDSS
     175      stars[N].average.R         = ra[i] + offsetRa[N] / 3600.0;
    162176      stars[N].average.D         = dec[i] + offsetDec[N] / 3600.0;
    163177      stars[N].average.dR        = NAN;
    … …  
    168182      stars[N].measure.dXccd     = ShortPixels(colcErr[N]);
    169183      stars[N].measure.dYccd     = ShortPixels(rowcErr[N]);
    170       // XXXYYY again the question what ZeroPt-zeropt[j] evalutes to - according to what it says above, this
     184      // XXXJester: again the question what ZeroPt-zeropt[j] evalutes to - according to what it says above, this
    171185      // should take .M back into instrumental magnitudes, but does it? In particular, the *actual* zeropoint
    172186      // that was used for the computation of this calibrated magnitude is in the tsField file; in general, it
    … …  
    176190      stars[N].measure.dM        = psfCountsErr[N];
    177191      stars[N].measure.Map       = fiberCounts[N] + ZeroPt - zeropt[j];
    178       stars[N].measure.Sky       = sky[N]; // adjust this to counts? XXXYYY It's given in maggies; to get back to counts,
     192      stars[N].measure.Sky       = sky[N]; // adjust this to counts? XXXJester: The SDSS Sky is given in maggies/arcsec^2; a maggie is a linear flux measure, 1 maggie corresponds to 0 magnitude, a surface brightness of 20 mag/square arcsec corresponds to 10-8 maggies per square arcsec.
    179193      stars[N].measure.dSky      = skyErr[N];
    180       stars[N].measure.FWx       = ShortPixels(seeing[j]); // reported in arcsec? XXXYYY: Yes!
     194      stars[N].measure.FWx       = ShortPixels(seeing[j]); // reported in arcsec? XXXJester:: Yes!
    181195      stars[N].measure.FWy       = ShortPixels(seeing[j]);
    182196      stars[N].measure.psfChisq  = prob_psf[N]; // XXX not really the correct value...
    183197      stars[N].measure.detID     = N;
    184198      stars[N].measure.t         = tzero[j] + clockRate*rowc[N]; // time since row 0
    185       stars[N].measure.dt        = 53.907456; // is this 2048*clockRate ? XXXYYY Nearly - 4 milliseconds more. shouldn't this be a #define?
     199      stars[N].measure.dt        = 53.907456; // is this 2048*clockRate ? XXXJester: Nearly - 4 milliseconds more. shouldn't this be a #define?
    186200
    187201      SetSDSSFlags (&stars[N], flags[N], flags2[N]);
    … …  
    199213      altaz (&alt, &az, sidtime - stars[N].average.R, stars[N].average.D, Latitude);
    200214
    201       // XXXYYY airmass is in tsField file
     215      // XXXJester: airmass is in tsField file
    202216      stars[N].measure.airmass   = 1.0 / dCOS(90.0 - alt);
    203217      stars[N].measure.az        = az;
    … …  
    206220  }
    207221
    208   // Throw away unneeded trailing (I hope) part of Stars
    209   REALLOCATE (stars, Stars, NFILTER*Nstars_used);
    210222
    211223  for (i = 0; i < NFILTER; i++) {
    … …  
    215227    // XXX for now, we define a totally fake coordinate system centered on the first listed star
    216228
    217     // XXXYYY what's happening here? *Exact* astrometry for SDSS can only be retrieved using the
     229    // XXXJester: what's happening here? *Exact* astrometry for SDSS can only be retrieved using the
    218230    // colour-dependent distortion coefficients from the tsField file. See
    219231    // www.sdss.org/dr7/products/general/astrometry.html
    … …  
    242254
    243255    images[N].NX = 2048;
    244     // XXXYYY NZ should be 1361, since the last 128 rows are the duplicate of the neighbouring frame, and only
     256    images[N].NY = 1490;
     257    // XXXJester: NY should be 1361, since the last 128 rows are the duplicate of the neighbouring frame, and only
    245258    // objects in 64 <= row < 1361+64 (or maybe the <= is at the upper limit?) are counted as belonging to
    246259    // "this" field.  I.e. we also need to read the status flag and use only primary and secondary objects
    … …  
    261274    // run. This is desired in the case of dvo. Hence we need exactly those objects with OK_RUN set, not more,
    262275    // not less.
    263 
    264     // And once more, NY = 1361, which led into the whole 'status' discussion. In any case, why is this 1490
    265     // not 1489?
    266     images[N].NY = 1490;
     276    //
     277    // And once more, I think we need NY = 1361, which led into the whole 'status' discussion. In any case,
     278    // why is this 1490 not 1489?
    267279
    268280    images[N].tzero = tzero[i];
    … …  
    270282
    271283    // set photcodes for the 5 images (SDSS_U,G,R,I,Z)
    272     // XXXYYY Why these? and not the 'u_sdss etc. 'ref' photcodes (see above where it says U_SDSS etc.)
     284    //
     285    // We use the 'dep' photcodes and not 'ref' because we can in principle recalibrate the SDSS photometry in
     286    // dvo; looks like 'ref' is reserved for real standard-star observations
    273287    images[N].photcode = photcode[i];
    274288
    … …  
    343357int SetSDSSFlags (Stars *star, unsigned int flags1, unsigned int flags2) {
    344358
    345   // XXX this is wrong, need to roll left to set the correct bit
    346359  if (flags1 & 0x00000002) star[0].measure.photFlags |= 0x0001; // BRIGHT            - 1  1
    347360  if (flags1 & 0x00000004) star[0].measure.photFlags |= 0x0002; // EDGE              - 1  2
Note: See TracChangeset for help on using the changeset viewer.