Index: /branches/eam_branches/ipp-20111122/Ohana/src/libdvo/include/dvo.h =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/libdvo/include/dvo.h (revision 33373) +++ /branches/eam_branches/ipp-20111122/Ohana/src/libdvo/include/dvo.h (revision 33374) @@ -142,5 +142,5 @@ */ -/* Average.code values -- these values are 32 bit (as of PS1_V1) */ +/* Average.flags values -- these values are 32 bit (as of PS1_V1) */ typedef enum { ID_STAR_FEW = 0x00000001, // used within relphot: skip star @@ -169,4 +169,15 @@ } DVOAverageFlags; +/* Secfilt.flags values -- these values are 32 bit (as of PS1_V1) */ +typedef enum { + ID_SECF_STAR_FEW = 0x00000001, // used within relphot: skip star + ID_SECF_STAR_POOR = 0x00000002, // used within relphot: skip star + ID_PHOTOM_PASS_0 = 0x00000100, // average magnitude calculated in 0th pass + ID_PHOTOM_PASS_1 = 0x00000200, // average magnitude calculated in 1th pass + ID_PHOTOM_PASS_2 = 0x00000400, // average magnitude calculated in 2th pass + ID_PHOTOM_PASS_3 = 0x00000800, // average magnitude calculated in 3th pass + ID_PHOTOM_PASS_4 = 0x00001000, // average magnitude calculated in 4th pass +} DVOSecfiltFlags; + /*** general dvo structures (internal use only / not IO) ***/ Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/include/relphot.h =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/include/relphot.h (revision 33373) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/include/relphot.h (revision 33374) @@ -216,4 +216,5 @@ float getMmos PROTO((off_t meas, int cat)); float getMrel PROTO((Catalog *catalog, off_t meas, int cat)); +short getUbercalDist PROTO((off_t meas, int cat)); Image *getimage PROTO((off_t N)); Image *getimages PROTO((off_t *N, off_t **LineNumber)); @@ -265,5 +266,5 @@ int setMrel PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); void setMrelFinal PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); -int setMrelOutput PROTO((Catalog *catalog, int Ncatalog, int mark, FlatCorrectionTable *flatcorr)); +int setMrelOutput PROTO((Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr)); int setMave PROTO((Catalog *catalog, int Ncatalog)); void set_ZP PROTO((double ZERO)); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/ImageOps.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/ImageOps.c (revision 33373) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/ImageOps.c (revision 33374) @@ -328,4 +328,18 @@ } +// returns image.Mcal - ff(x,y) +short getUbercalDist (off_t meas, int cat) { + + off_t i; + short distance; + + i = MeasureToImage[cat][meas]; + if (i == -1) return (1000); + + if (image[i].flags & IMAGE_BAD) return (1000); + distance = image[i].dummy3; // XXX dummy3 == ubercalDist + return (distance); +} + float getMflat (off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog) { @@ -393,6 +407,12 @@ } - /* we optionally do not recalibrate images with UBERCAL zero points */ - if (KEEP_UBERCAL && (image[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; + // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0: + // we optionally do not recalibrate images with UBERCAL zero points + if (image[i].flags & ID_IMAGE_PHOTOM_UBERCAL) { + image[i].dummy3 = 0; // dummy3 = ubercalDist + if (KEEP_UBERCAL) continue; + } + + int minUbercalDist = 1000; N = 0; @@ -436,4 +456,6 @@ } + minUbercalDist = MIN (catalog[c].secfilt[n*Nsecfilt].M_80, minUbercalDist); // M_80 == ubercalDist + // fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f %6.3f %6.3f %6.3f %6.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, Msys, Mrel, Mmos, Mgrid, Mflat); @@ -473,4 +495,6 @@ image[i].Mcal = 0.0; } + + image[i].dummy3 = minUbercalDist; } free (list); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/StarOps.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/StarOps.c (revision 33373) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/StarOps.c (revision 33374) @@ -61,4 +61,6 @@ // XXX accumulate all secfilt values in a single pass? + + int minUbercalDist = 1000; int Ns; @@ -102,4 +104,7 @@ list[N] = Msys - Mcal - Mmos - Mgrid; + int myUbercalDist = getUbercalDist(m,i); + minUbercalDist = MIN(minUbercalDist, myUbercalDist); + // dlist gives the error, which is used as the weight in WT_MEAN. // we can modify the resulting weight in a few ways: @@ -137,8 +142,10 @@ liststats (list, dlist, N, &stats); - + catalog[i].secfilt[Nsecfilt*j+Nsec].M = stats.mean; catalog[i].secfilt[Nsecfilt*j+Nsec].dM = stats.sigma; catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq + 1e-4) : NAN_S_SHORT; + + catalog[i].secfilt[Nsecfilt*j+Nsec].M_80 = minUbercalDist; // M_80 == ubercalDist } } @@ -159,5 +166,5 @@ // * setMave updates 2MASS average flags // * setMave updates average EXT flags (PS1 and 2MASS) -int setMrelOutput (Catalog *catalog, int Ncatalog, int mark, FlatCorrectionTable *flatcorr) { +int setMrelOutput (Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr) { int i, N; @@ -192,4 +199,6 @@ if (catalog[i].found[Nsecfilt*j+Nsec]) continue; + int Ncode = 0; + N = 0; m = catalog[i].averageT[j].measureOffset; @@ -200,10 +209,8 @@ if (!code) continue; if (code->equiv != thisCode) { continue; } + Ncode ++; if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; - // XXX allow REF stars (no Image Entry) to be included in the calculation this - // should be optionally set, and should allow for REF stars to be downweighted by - // more than their reported errors. how such info is carried is unclear... if (getImageEntry (m, i) < 0) { Mcal = Mmos = Mgrid = 0; @@ -247,6 +254,25 @@ if (N < 1) continue; + // force WT_MEAN or MEAN here? + liststats (list, dlist, N, &stats); - if (mark) catalog[i].found[Nsecfilt*j+Nsec] = TRUE; + catalog[i].found[Nsecfilt*j+Nsec] = TRUE; + switch (pass) { + case 0: + catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_0; + break; + case 1: + catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_1; + break; + case 2: + catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_2; + break; + case 3: + catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_3; + break; + case 4: + catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_4; + break; + } /* use sigma or error in dM for output? */ @@ -254,4 +280,6 @@ catalog[i].secfilt[Nsecfilt*j+Nsec].dM = MAX (stats.error, stats.sigma); catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT; + catalog[i].secfilt[Nsecfilt*j+Nsec].Ncode = Ncode; + catalog[i].secfilt[Nsecfilt*j+Nsec].Nused = stats.Nmeas; } } @@ -540,10 +568,10 @@ for (Ns = 0; Ns < Nphotcodes; Ns++) { - /* on final processing, skip stars already measured */ - if (final && catalog[i].found[j]) continue; - int thisCode = photcodes[Ns][0].code; int Nsec = GetPhotcodeNsec(thisCode); + /* on final processing, skip stars already measured */ + if (final && catalog[i].found[Nsecfilt*j + Nsec]) continue; + /* skip bad stars to prevent them from becoming good (on inner sample) */ if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; @@ -557,6 +585,5 @@ if (ecode != thisCode) { continue; } - // why am I skipping this?? - /* if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; */ + // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers Mcal = getMcal (m, i, flatcorr, catalog); if (isnan(Mcal)) { Ncal ++; continue; } @@ -602,5 +629,5 @@ if (ecode != thisCode) { continue; } - /* if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; */ + // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers Mcal = getMcal (m, i, flatcorr, catalog); if (isnan(Mcal)) continue; @@ -725,7 +752,4 @@ list[n] = pow (10.0, 0.01*Xm); dlist[n] = 1; - if (list[n] > 100000) { - fprintf (stderr, "!"); - } n++; } @@ -765,7 +789,4 @@ dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dM; if (isnan(dM)) continue; - if (dM > 100) { - fprintf (stderr, "."); - } list[n] = dM; dlist[n] = 1; Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/select_images.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/select_images.c (revision 33373) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/select_images.c (revision 33374) @@ -239,9 +239,11 @@ } if (RESET) { - // XXX assignMcal (&image[nimage], (double *) NULL, -1); // XXX this needs to be thought through a bit more - image[nimage].Mcal = 0.0; - image[nimage].dMcal = NAN; + if (!KEEP_UBERCAL) { + image[nimage].Mcal = 0.0; + image[nimage].dMcal = NAN; + } image[nimage].flags &= ~ID_IMAGE_PHOTOM_POOR; + image[nimage].dummy3 = 1000; // dummy3 == ubercalDist } line_number[nimage] = i; Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/setMrelFinal.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/setMrelFinal.c (revision 33373) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/setMrelFinal.c (revision 33374) @@ -21,4 +21,12 @@ ID_MEAS_AREA; // detetion was outside acceptable area of device + // flags used by the photometry analysis (excluding UBERCAL) + unsigned int PHOTOM_PASSES = + ID_PHOTOM_PASS_0 | // average measured at pass 0 + ID_PHOTOM_PASS_1 | // average measured at pass 1 + ID_PHOTOM_PASS_2 | // average measured at pass 2 + ID_PHOTOM_PASS_3 | // average measured at pass 3 + ID_PHOTOM_PASS_4; // average measured at pass 4 + int Ns; for (Ns = 0; Ns < Nphotcodes; Ns++) { @@ -31,4 +39,8 @@ catalog[0].secfilt[Nsecfilt*i+Nsec].dM = NAN; catalog[0].secfilt[Nsecfilt*i+Nsec].Xm = NAN_S_SHORT; + catalog[0].secfilt[Nsecfilt*j+Nsec].Ncode = 0; + catalog[0].secfilt[Nsecfilt*j+Nsec].Nused = 0; + + catalog[0].secfilt[Nsecfilt*i+Nsec].flags &= ~PHOTOM_PASSES; m = catalog[0].average[i].measureOffset; @@ -64,16 +76,14 @@ } + clean_measures (catalog, 1, TRUE, flatcorr); /* mark outliers ID_MEAS_POOR_PHOTOM */ for (i = 0; i < 5; i++) { skip_measurements (catalog, i, flatcorr); /* set ID_MEAS_SKIP for measures to be skipped */ - setMrelOutput (catalog, 1, FALSE, flatcorr); /* set Mrel using all measures not skipped */ - clean_measures (catalog, 1, TRUE, flatcorr); /* mark outliers ID_MEAS_POOR_PHOTOM */ - setMrelOutput (catalog, 1, TRUE, flatcorr); /* set Mrel using remaining measures */ + setMrelOutput (catalog, 1, i, flatcorr); /* set Mrel using remaining measures */ } setMcalOutput (catalog, 1, flatcorr); - - int Nsecfilt = GetPhotcodeNsecfilt (); /* clear ID_STAR_POOR, ID_STAR_FEW values before writing ??? */ /* ID_MEAS_NOCAL is an internal bit, so it should be cleared */ + // XXX eventually: move these internal flags to an internal-only catalog->FOO element (like found) for (i = 0; i < catalog[0].Naverage; i++) { catalog[0].average[i].flags &= ~ID_STAR_FEW; @@ -95,7 +105,8 @@ // acceptance on each pass. The rejections at the different passes are: -// 0 : only use measurements thought to be not extended and dM < SIGMA_LIM (REMOVED) -// 1 : use extended and low-sig detections (REMOVED) -// 3 : accept the measurements marked as outliers; accept images marked as outliers +// 0 : only use measurements thought to be GOOD (photflags not POOR or BAD) +// 1 : accept measurements thought to be POOR based on photflags +// 2 : accept the measurements marked as outliers; accept images marked as outliers +// 3 : accept measurements thought to be BAD based on photflags // 4 : accept the measurements outside of the instrumental magnitude limits (eg, SAT) @@ -106,83 +117,88 @@ off_t i, k, m; - int ecode, found, Ns; off_t Ntot, Ntry, Nkeep, Nskip; float mag; + int Nsecfilt = GetPhotcodeNsecfilt (); + Ntot = Ntry = Nskip = Nkeep = 0; /* allow measures from images marked POOR and FEW */ - if (pass >= 3) IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL; + if (pass >= 2) IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL; /* allow measures marked as outliers (POOR) and off image region (AREA) */ - if (pass >= 3) MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM; - - /* mark measures which should be ignored on second pass */ + if (pass >= 2) MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM; + + /* mark measures which should be ignored on this pass */ for (i = 0; i < catalog[0].Naverage; i++) { Ntot += catalog[0].average[i].Nmeasure; - // this function does not check secfilt entries separately, but there is no harm in - // twiddling the so skip a star here only we can only - // have we gotten - - if (catalog[0].found[i]) continue; - XXX : this should be done by secfilt. think about this... - - m = catalog[0].average[i].measureOffset; - for (k = 0; k < catalog[0].average[i].Nmeasure; k++, m++) { - Ntry++; - - /* clear SKIP for all measures at first */ - catalog[0].measureT[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; - - /** never use these measurements (wrong photcode, bad time range) */ - /* skipped via NOCAL, don't mark as skipped */ - ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].photcode); - found = FALSE; - for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) { - if (ecode == photcodes[Ns][0].code) found = TRUE; + // mark measurements for each secfilt separately + int Ns; + for (Ns = 0; Ns < Nphotcodes; Ns++) { + + int thisCode = photcodes[Ns][0].code; + int Nsec = GetPhotcodeNsec(thisCode); + + /* star/photcodes already calibrated */ + if (catalog[0].found[Nsecfilt*i+Nsec]) continue; + + m = catalog[0].average[i].measureOffset; + for (k = 0; k < catalog[0].average[i].Nmeasure; k++, m++) { + Ntry++; + + // skip measurements not related to this photcode + PhotCode *code = GetPhotcodebyCode (catalog[0].measure[m].photcode); + if (!code) continue; + if (code->equiv != thisCode) continue; + + // skip measurements by time range (mark as skipped for this secfilt) + // XXX note that this is a bit dangerous : some objects may never get calibrated + if (TimeSelect) { + if (catalog[0].measure[m].t < TSTART) goto skip; + if (catalog[0].measure[m].t > TSTOP) goto skip; + } + + // we now have a measurement relavant to this photcode (and time range) + + // clear SKIP for all measures at first + catalog[0].measureT[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; + + // skip measurements from BAD images and mosaics + // do NOT skip measurements without a matching image + if (isnan(getMcal (m, 0, flatcorr, catalog))) goto skip; + if (isnan(getMmos (m, 0))) goto skip; + + // PASS 4 : skip measurements by inst mag limit + if ((pass < 4) && ImagSelect) { + mag = PhotInst (&catalog[0].measure[m]); + if (mag < ImagMin) goto skip; + if (mag > ImagMax) goto skip; + } + + // PASS 3 : accept bad measurements (eg, SAT, CR) + if (pass < 3) { + if (catalog[0].measure[m].photFlags & code->photomBadMask) goto skip; + if ((catalog[i].measure[m].photcode > 10000) && (catalog[i].measure[m].photcode < 10500)) { + if (catalog[i].measure[m].psfQual < 0.85) goto skip; + } + } + + // PASS 2 : internal outliers accepted + + // PASS 1 : accept poor measurements as well (eg, POOR FIT, etc) + if ((pass < 1) && (catalog[0].measure[m].photFlags & code->photomPoorMask)) goto skip; + + // PASS 0 : skip poor, outlier, bad measurements as well (eg, POOR FIT, etc) + + Nkeep ++; + + continue; + + skip: + catalog[0].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM; + catalog[0].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM; + Nskip ++; } - if (!found) continue; - - /* skip measurements by time range */ - if (TimeSelect) { - if (catalog[0].measure[m].t < TSTART) goto skip; - if (catalog[0].measure[m].t > TSTOP) goto skip; - } - - /* skip measurements with sat. dophot values */ - // XXX we need to fix this: use the photFlags to ID saturation... - // if ((pass < 4) && (catalog[0].measure[m].dophot == 10)) goto skip; - - /* skip measurements from BAD images and mosaics */ - /* do NOT skip measurements without a matching image */ - if (isnan(getMcal (m, 0, flatcorr, catalog))) goto skip; - if (isnan(getMmos (m, 0))) goto skip; - - /* skip measurements by inst mag limit */ - if ((pass < 4) && ImagSelect) { - mag = PhotInst (&catalog[0].measure[m]); - if (mag < ImagMin) goto skip; - if (mag > ImagMax) goto skip; - } - - // d1 = (catalog[0].measure[m].dophot == 1); - // d2 = (catalog[0].measure[m].dophot == 2); - // fprintf (stderr, "fix the dophot / type selections!"); - // abort (); - - /* skip measurements by measurement error */ - // if ((pass < 2) && !(d1 || d2)) goto skip; - - /* if ((pass < 1) && (catalog[0].measure[m].dM > SIGMA_LIM)) goto skip; */ - // if ((pass < 1) && !d1) goto skip; - Nkeep ++; - - continue; - - skip: - catalog[0].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM; - catalog[0].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM; - Nskip ++; } }