IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 24980


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

Comments (marked by XXXJester) about possibly reporting both RMS and sigma, rather than MAX(stats.error,stats.sigma)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/sj_branches/sj_SDSSaddstar_branch_r23928_20090420/Ohana/src/relphot/src/StarOps.c

    r21508 r24980  
    88
    99  int i, j;
    10  
     10
    1111  Nmax = 0;
    1212  for (i = 0; i < Ncatalog; i++) {
    … …  
    1818  ALLOCATE (list, double, MAX (1, Nmax));
    1919  ALLOCATE (dlist, double, MAX (1, Nmax));
    20 } 
     20}
    2121
    2222float getMrel (Catalog *catalog, int meas, int cat) {
    … …  
    2626
    2727  ave = catalog[cat].measure[meas].averef;
    28   if (catalog[cat].average[ave].flags & STAR_BAD) return (NAN); 
    29  
     28  if (catalog[cat].average[ave].flags & STAR_BAD) return (NAN);
     29
    3030  value = catalog[cat].secfilt[PhotNsec*ave+PhotSec].M;
    3131  return (value);
    … …  
    4444
    4545      /* calculate the average value for a single star */
    46       if (catalog[i].average[j].flags & STAR_BAD) continue; 
     46      if (catalog[i].average[j].flags & STAR_BAD) continue;
    4747      m = catalog[i].average[j].measureOffset;
    4848
    4949      N = 0;
    5050      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
    51         if (catalog[i].measure[m].dbFlags & MEAS_BAD) {
    52           Nbad ++;
    53           continue;
    54         }
    55         // XXX allow REF stars (no Image Entry) to be included in the calculation this
    56         // should be optionally set, and should allow for REF stars to be downweighted by
    57         // more than their reported errors.  how such info is carried is unclear...
    58         if (getImageEntry (m, i) < 0) {
    59           Mcal = Mmos = Mgrid = 0;
    60         } else {
    61           Mcal  = getMcal  (m, i);
    62           if (isnan(Mcal)) {
    63             Ncal ++;
    64             continue;
    65           }
    66           Mmos  = getMmos  (m, i);
    67           if (isnan(Mmos)) {
    68             Nmos ++;
    69             continue;
    70           }
    71           Mgrid = getMgrid (m, i);
    72           if (isnan(Mgrid)) {
    73             Ngrid++;
    74             continue;
    75           }
    76         }
    77 
    78         Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
    79         if (isnan(Msys)) {
    80           Nsys++;
    81           continue;
    82         }
    83         list[N] = Msys - Mcal - Mmos - Mgrid;
    84         dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
    85         N++;
     51        if (catalog[i].measure[m].dbFlags & MEAS_BAD) {
     52          Nbad ++;
     53          continue;
     54        }
     55        // XXX allow REF stars (no Image Entry) to be included in the calculation this
     56        // should be optionally set, and should allow for REF stars to be downweighted by
     57        // more than their reported errors.  how such info is carried is unclear...
     58        if (getImageEntry (m, i) < 0) {
     59          Mcal = Mmos = Mgrid = 0;
     60        } else {
     61          Mcal  = getMcal  (m, i);
     62          if (isnan(Mcal)) {
     63            Ncal ++;
     64            continue;
     65          }
     66          Mmos  = getMmos  (m, i);
     67          if (isnan(Mmos)) {
     68            Nmos ++;
     69            continue;
     70          }
     71          Mgrid = getMgrid (m, i);
     72          if (isnan(Mgrid)) {
     73            Ngrid++;
     74            continue;
     75          }
     76        }
     77
     78        Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
     79        if (isnan(Msys)) {
     80          Nsys++;
     81          continue;
     82        }
     83        list[N] = Msys - Mcal - Mmos - Mgrid;
     84        dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
     85        N++;
    8686      }
    8787
    … …  
    8989
    9090      if (N <= STAR_TOOFEW) { /* too few measurements */
    91         catalog[i].average[j].flags |= ID_STAR_FEW;
    92         Nfew ++;
     91        catalog[i].average[j].flags |= ID_STAR_FEW;
     92        Nfew ++;
    9393      } else {
    94         catalog[i].average[j].flags &= ~ID_STAR_FEW;
    95       } 
     94        catalog[i].average[j].flags &= ~ID_STAR_FEW;
     95      }
    9696
    9797      liststats (list, dlist, N, &stats);
    … …  
    128128
    129129      /* skip stars already calibrated */
    130       if (catalog[i].found[j]) continue; 
     130      if (catalog[i].found[j]) continue;
    131131
    132132      N = 0;
    133133      m = catalog[i].average[j].measureOffset;
    134134      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
    135         if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
    136         // XXX allow REF stars (no Image Entry) to be included in the calculation this
    137         // should be optionally set, and should allow for REF stars to be downweighted by
    138         // more than their reported errors.  how such info is carried is unclear...
    139         if (getImageEntry (m, i) < 0) {
    140           Mcal = Mmos = Mgrid = 0;
    141         } else {
    142           Mcal  = getMcal  (m, i);
    143           if (isnan(Mcal)) continue;
    144           Mmos  = getMmos  (m, i);
    145           if (isnan(Mmos)) continue;
    146           Mgrid = getMgrid (m, i);
    147           if (isnan(Mgrid)) continue;
    148         }
    149 
    150         Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
    151         list[N] = Msys - Mcal - Mmos - Mgrid;
    152         dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
    153         N++;
     135        if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
     136        // XXX allow REF stars (no Image Entry) to be included in the calculation this
     137        // should be optionally set, and should allow for REF stars to be downweighted by
     138        // more than their reported errors.  how such info is carried is unclear...
     139        if (getImageEntry (m, i) < 0) {
     140          Mcal = Mmos = Mgrid = 0;
     141        } else {
     142          Mcal  = getMcal  (m, i);
     143          if (isnan(Mcal)) continue;
     144          Mmos  = getMmos  (m, i);
     145          if (isnan(Mmos)) continue;
     146          Mgrid = getMgrid (m, i);
     147          if (isnan(Mgrid)) continue;
     148        }
     149
     150        Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
     151        list[N] = Msys - Mcal - Mmos - Mgrid;
     152        dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
     153        N++;
    154154      }
    155155      if (N < 1) continue;
    … …  
    195195      for (Ns = 0; Ns < Nsecfilt; Ns++) {
    196196
    197         code = GetPhotcodebyNsec (Ns);
    198         Nc = code[0].code;
    199        
    200         N = 0;
    201         m = catalog[i].average[j].measureOffset;
    202         for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
    203           if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
    204           if (GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode) != Nc) continue;
    205 
    206           Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
    207           if (isnan(Msys)) continue;
    208 
    209           // XXX this is a hack for the 2MASS search; better to save an average value?
    210           if (catalog[i].measure[m].psfQual < 0.85) continue;
    211 
    212           list[N] = Msys - catalog[i].measure[m].Mcal;
    213           dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
    214           N++;
    215         }
    216         if (N < 1) continue;
    217 
    218         liststats (list, dlist, N, &stats);
    219 
    220         /* use sigma or error in dM for output? */
    221         catalog[i].secfilt[Nsecfilt*j+Ns].M  = stats.mean;
    222         catalog[i].secfilt[Nsecfilt*j+Ns].dM = MAX (stats.error, stats.sigma);
    223         catalog[i].secfilt[Nsecfilt*j+Ns].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
    224         catalog[i].secfilt[Nsecfilt*j+Ns].Ncode = N;
    225         catalog[i].secfilt[Nsecfilt*j+Ns].Nused = stats.Nmeas;
     197        code = GetPhotcodebyNsec (Ns);
     198        Nc = code[0].code;
     199
     200        N = 0;
     201        m = catalog[i].average[j].measureOffset;
     202        for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
     203          if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
     204          if (GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode) != Nc) continue;
     205
     206          Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
     207          if (isnan(Msys)) continue;
     208
     209          // XXX this is a hack for the 2MASS search; better to save an average value?
     210          if (catalog[i].measure[m].psfQual < 0.85) continue;
     211
     212          list[N] = Msys - catalog[i].measure[m].Mcal;
     213          dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
     214          N++;
     215        }
     216        if (N < 1) continue;
     217
     218        liststats (list, dlist, N, &stats);
     219
     220        /* use sigma or error in dM for output? */
     221        // XXXJester: why not report *both* the RMS and the sigma, where one is the error-propagated standard
     222        // error on the mean, and the other is the RMS of the values about the mean [with or without scaling
     223        // by sqrt(N)]? If the two diverge by a lot, we have some un-understood source of systematic errors
     224        // that could be QA-ed by looking for this divergence, but will otherwise only be found if people
     225        // think actively about the expectations for what the errors should be, in effect doing this
     226        // calculation by hand.
     227        catalog[i].secfilt[Nsecfilt*j+Ns].M  = stats.mean;
     228        // XXXJester: Neither of stats.error,stats.sigma must be NaN for the MAX() comparison to be meaningful!
     229        catalog[i].secfilt[Nsecfilt*j+Ns].dM = MAX (stats.error, stats.sigma);
     230        catalog[i].secfilt[Nsecfilt*j+Ns].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
     231        catalog[i].secfilt[Nsecfilt*j+Ns].Ncode = N;
     232        catalog[i].secfilt[Nsecfilt*j+Ns].Nused = stats.Nmeas;
    226233      }
    227234    }
    … …  
    247254      m = catalog[i].average[j].measureOffset;
    248255      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
    249         if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
    250         Mcal  = getMcal  (m, i);
    251         if (isnan(Mcal)) continue;
    252         Mmos  = getMmos  (m, i);
    253         if (isnan(Mmos)) continue;
    254         Mgrid = getMgrid (m, i);
    255         if (isnan(Mgrid)) continue;
    256         catalog[i].measure[m].Mcal = Mcal + Mmos + Mgrid;
     256        if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
     257        Mcal  = getMcal  (m, i);
     258        if (isnan(Mcal)) continue;
     259        Mmos  = getMmos  (m, i);
     260        if (isnan(Mmos)) continue;
     261        Mgrid = getMgrid (m, i);
     262        if (isnan(Mgrid)) continue;
     263        catalog[i].measure[m].Mcal = Mcal + Mmos + Mgrid;
    257264      }
    258265    }
    … …  
    273280  /* find Xm median -> ChiSq lim must be > median */
    274281  for (i = Ntot = 0; i < Ncatalog; i++) {
    275     Ntot += catalog[i].Naverage; 
     282    Ntot += catalog[i].Naverage;
    276283  }
    277284  ALLOCATE (xlist, double, Ntot);
    … …  
    290297    }
    291298  }
    292  
     299
    293300  initstats ("MEAN");
    294301  liststats (xlist, dlist, Ntot, &stats);
    … …  
    306313      mark = (dM > MaxScatter) || (Xm == NAN_S_SHORT) || (Chisq > MaxChisq);
    307314      if (mark) {
    308         catalog[i].average[j].flags |= ID_STAR_POOR;
    309         Ndel ++;
     315        catalog[i].average[j].flags |= ID_STAR_POOR;
     316        Ndel ++;
    310317      } else {
    311         catalog[i].average[j].flags &= ~ID_STAR_POOR;
     318        catalog[i].average[j].flags &= ~ID_STAR_POOR;
    312319      }
    313320      Nave ++;
    … …  
    341348  ALLOCATE (ilist, int, Nmax);
    342349  ALLOCATE (tlist, double, Nmax);
    343  
     350
    344351  /* it makes no sense to mark 3-sigma outliers with <5 measurements */
    345352  TOOFEW = MAX (5, STAR_TOOFEW);
    … …  
    351358
    352359      /* skip bad stars to prevent them from becoming good (on inner sample) */
    353       if (catalog[i].average[j].flags & STAR_BAD) continue; 
     360      if (catalog[i].average[j].flags & STAR_BAD) continue;
    354361
    355362      /* on final processing, skip stars already measured */
    356       if (final && catalog[i].found[j]) continue; 
     363      if (final && catalog[i].found[j]) continue;
    357364
    358365      /* accumulate list of valid measurements */
    … …  
    360367      N = 0;
    361368      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
    362         /* if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue; */
    363         Mcal  = getMcal  (m, i);
    364         if (isnan(Mcal)) continue;
    365         Mmos  = getMmos  (m, i);
    366         if (isnan(Mmos)) continue;
    367         Mgrid = getMgrid (m, i);
    368         if (isnan(Mgrid)) continue;
    369 
    370         Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
    371         list[N] = Msys - Mcal - Mmos - Mgrid;
    372         dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
    373         N++;
     369        /* if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue; */
     370        Mcal  = getMcal  (m, i);
     371        if (isnan(Mcal)) continue;
     372        Mmos  = getMmos  (m, i);
     373        if (isnan(Mmos)) continue;
     374        Mgrid = getMgrid (m, i);
     375        if (isnan(Mgrid)) continue;
     376
     377        Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
     378        list[N] = Msys - Mcal - Mmos - Mgrid;
     379        dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
     380        N++;
    374381      }
    375382      if (N <= TOOFEW) continue;
    … …  
    380387      stats.sigma = MAX (MIN_ERROR, stats.sigma); /* if measurements agree too well, sigma -> 0.0 */
    381388      for (k = m = 0; k < N; k++) {
    382         if (fabs (list[k] - stats.median) < Ns*stats.sigma) {
    383           list[m] = list[k];
    384           m++;
    385         }
     389        if (fabs (list[k] - stats.median) < Ns*stats.sigma) {
     390          list[m] = list[k];
     391          m++;
     392        }
    386393      }
    387394      initstats ("MEAN");
    … …  
    395402      N = 0;
    396403      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
    397         /* if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue; */
    398         Mcal  = getMcal  (m, i);
    399         if (isnan(Mcal)) continue;
    400         Mmos  = getMmos  (m, i);
    401         if (isnan(Mmos)) continue;
    402         Mgrid = getMgrid (m, i);
    403         if (isnan(Mgrid)) continue;
    404 
    405         Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
    406         list[N] = Msys - Mcal - Mmos - Mgrid;
    407         dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
    408         ilist[N] = m;
    409         N++;
    410         Nave ++;
     404        /* if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue; */
     405        Mcal  = getMcal  (m, i);
     406        if (isnan(Mcal)) continue;
     407        Mmos  = getMmos  (m, i);
     408        if (isnan(Mmos)) continue;
     409        Mgrid = getMgrid (m, i);
     410        if (isnan(Mgrid)) continue;
     411
     412        Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
     413        list[N] = Msys - Mcal - Mmos - Mgrid;
     414        dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
     415        ilist[N] = m;
     416        N++;
     417        Nave ++;
    411418      }
    412419      if (N < TOOFEW) continue;
    … …  
    414421      /* mark bad measures */
    415422      for (k = 0; k < N; k++) {
    416         if (fabs (list[k] - stats.median) > Ns*stats.sigma) {
    417           catalog[i].measure[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM;
    418           Ndel ++;
    419         }
     423        if (fabs (list[k] - stats.median) > Ns*stats.sigma) {
     424          catalog[i].measure[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM;
     425          Ndel ++;
     426        }
    420427      }
    421428      IMAGE_BAD = image_bad;
    … …  
    437444  Ntot = 0;
    438445  for (i = 0; i < Ncatalog; i++) {
    439     Ntot += catalog[i].Naverage; 
     446    Ntot += catalog[i].Naverage;
    440447  }
    441448
    … …  
    448455
    449456      /* calculate the average value for a single star */
    450       if (catalog[i].average[j].flags & STAR_BAD) continue; 
     457      if (catalog[i].average[j].flags & STAR_BAD) continue;
    451458      m = catalog[i].average[j].measureOffset;
    452459
    453460      N = 0;
    454461      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
    455         Mcal = getMcal  (m, i);
    456         if (isnan(Mcal)) continue;
    457         Mmos = getMmos  (m, i);
    458         if (isnan(Mmos)) continue;
    459         Mgrid = getMgrid (m, i);
    460         if (isnan(Mgrid)) continue;
    461         N++;
    462       }
    463      
     462        Mcal = getMcal  (m, i);
     463        if (isnan(Mcal)) continue;
     464        Mmos = getMmos  (m, i);
     465        if (isnan(Mmos)) continue;
     466        Mgrid = getMgrid (m, i);
     467        if (isnan(Mgrid)) continue;
     468        N++;
     469      }
     470
    464471      list[n] = N;
    465472      dlist[n] = 1;
    … …  
    482489  Ntot = 0;
    483490  for (i = 0; i < Ncatalog; i++) {
    484     Ntot += catalog[i].Naverage; 
     491    Ntot += catalog[i].Naverage;
    485492  }
    486493
    … …  
    493500
    494501      /* calculate the average value for a single star */
    495       if (catalog[i].average[j].flags & STAR_BAD) continue; 
     502      if (catalog[i].average[j].flags & STAR_BAD) continue;
    496503
    497504      Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
    … …  
    518525  Ntot = 0;
    519526  for (i = 0; i < Ncatalog; i++) {
    520     Ntot += catalog[i].Naverage; 
     527    Ntot += catalog[i].Naverage;
    521528  }
    522529
    … …  
    529536
    530537      /* calculate the average value for a single star */
    531       if (catalog[i].average[j].flags & STAR_BAD) continue; 
     538      if (catalog[i].average[j].flags & STAR_BAD) continue;
    532539
    533540      dM = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
    … …  
    559566  for (i = 0; i < Ncatalog; i++) {
    560567    for (j = 0; j < catalog[i].Naverage; j++) {
    561       if (catalog[i].average[j].flags & STAR_BAD) continue; 
     568      if (catalog[i].average[j].flags & STAR_BAD) continue;
    562569      dMrel = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
    563570      bin = dMrel / 0.00025;
    … …  
    613620  Graphdata graphdata;
    614621
    615   N = 0; 
     622  N = 0;
    616623  for (i = 0; i < Ncatalog; i++) {
    617624    N += catalog[i].Naverage;
Note: See TracChangeset for help on using the changeset viewer.