IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36524


Ignore:
Timestamp:
Feb 18, 2014, 2:27:35 PM (12 years ago)
Author:
eugene
Message:

independent lists for averaging psf, aper, kron mags

Location:
branches/eam_branches/ipp-20140214/Ohana/src/relphot
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r36491 r36524  
    8686  // others are only used on the final output steps.
    8787
    88   double *list;               // list of measure.mag values for a given star
    89   double *dlist;              // mag errors for a star
    90   double *wlist;              // weights to use for mean mags
    91   double *aplist;             // ap mags for a star
    92   double *kronlist;           // kron mags for a star
     88  double *Mpsflist;           // list of measure.mag values for a given star
     89  double *dpsflist;           // mag errors for a star
     90  double *wpsflist;           // weights to use for mean mags
     91
     92  double *Maplist;            // ap mags for a star
     93  double *daplist;            // ap mags for a star
     94  double *waplist;            // ap mags for a star
     95
     96  double *Mkronlist;          // kron mags for a star
    9397  double *dkronlist;          // kron mag errors
     98  double *wkronlist;          // kron mag errors
    9499
    95100  double *psfqf_list;         // psfqf for all filters
    … …  
    97102  double *stargal_list;       // stargal for all filters
    98103
    99   StatType stats;
     104  StatType psfstats;
    100105  StatType apstats;
    101106  StatType kronstats;

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

    r36491 r36524  
    6363  results->Ngrid = 0;
    6464  if (allocLists) {
    65     ALLOCATE (results->list, double, Nmax);
    66     ALLOCATE (results->dlist, double, Nmax);
    67     ALLOCATE (results->wlist, double, Nmax);
     65    ALLOCATE (results->Mpsflist, double, Nmax);
     66    ALLOCATE (results->dpsflist, double, Nmax);
     67    ALLOCATE (results->wpsflist, double, Nmax);
    6868  }
    6969}
    7070
    7171void SetMrelInfoFree (SetMrelInfo *results) {
    72   free (results->list);
    73   free (results->dlist);
    74   free (results->wlist);
     72  free (results->Mpsflist);
     73  free (results->dpsflist);
     74  free (results->wpsflist);
    7575}
    7676
    7777void SetMrelInfoAccum (SetMrelInfo *summary, SetMrelInfo *results) {
    7878  summary->Nfew  += results->Nfew ;
    79   summary->Ncode  += results->Ncode ;
     79  summary->Ncode += results->Ncode ;
    8080  summary->Nsys  += results->Nsys ;
    8181  summary->Nbad  += results->Nbad ;
    … …  
    158158  SetMrelInfoInit (&summary, FALSE);
    159159  SetMrelInfoInit (&results, TRUE); // allocates results->list,dlist,wlist
    160   ALLOCATE (results.aplist, double, Nmax);
    161   ALLOCATE (results.kronlist, double, Nmax);
     160
     161  ALLOCATE (results.Maplist,  double, Nmax);
     162  ALLOCATE (results.daplist,  double, Nmax);
     163  ALLOCATE (results.waplist,  double, Nmax);
     164
     165  ALLOCATE (results.Mkronlist, double, Nmax);
    162166  ALLOCATE (results.dkronlist, double, Nmax);
     167  ALLOCATE (results.wkronlist, double, Nmax);
    163168
    164169  ALLOCATE (results.psfqf_list, double, Nmax);
    … …  
    183188
    184189  SetMrelInfoFree (&results);
    185   free (results.aplist);
    186   free (results.kronlist);
     190  free (results.Maplist);
     191  free (results.daplist);
     192  free (results.waplist);
     193
     194  free (results.Mkronlist);
    187195  free (results.dkronlist);
     196  free (results.wkronlist);
    188197
    189198  free (results.psfqf_list);
    … …  
    331340  exit (3);
    332341
     342# if (0)
     343
    333344  off_t j, k, m, ID;
    334345  int N;
    … …  
    336347
    337348  StatType stats, apstats, kronstats;
    338   liststats_setmode (&stats, STATMODE);
     349  liststats_setmode (&psfstats, STATMODE);
    339350  liststats_setmode (&apstats, STATMODE);
    340351  liststats_setmode (&kronstats, STATMODE);
    … …  
    837848  if (primaryCell) free (primaryCell);
    838849  return (TRUE);
     850# endif
    839851}
    840852

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

    r36491 r36524  
    3434  off_t j;
    3535
    36   liststats_setmode (&results->stats, STATMODE);
     36  liststats_setmode (&results->psfstats, STATMODE);
    3737  liststats_setmode (&results->apstats, STATMODE);
    3838  liststats_setmode (&results->kronstats, STATMODE);
    … …  
    6767
    6868  off_t k;
    69   int N;
    7069  float Msys = 0, Mcal= 0, Mmos = 0, Mgrid = 0;
    7170
    72   double *list      = results->list;
    73   double *dlist     = results->dlist;
    74   double *wlist     = results->wlist;
    75   double *aplist    = results->aplist;
    76   double *kronlist  = results->kronlist;
     71  // we are measuring means for 3 types of mags: psf, ap, kron.  I am using the psf mag
     72  // error for the ap mags, but krons have their own errors.  it is an open question if I
     73  // should be doing weighted or unweighted fits (this is a user option)
     74  double *Mpsflist  = results->Mpsflist;
     75  double *dpsflist  = results->dpsflist;
     76  double *wpsflist  = results->wpsflist;
     77
     78  double *Maplist   = results->Maplist;
     79  double *daplist   = results->daplist;
     80  double *waplist   = results->waplist;
     81
     82  double *Mkronlist = results->Mkronlist;
    7783  double *dkronlist = results->dkronlist;
    78 
    79   StatType *stats = &results->stats;
    80   StatType *apstats = &results->apstats;
     84  double *wkronlist = results->wkronlist;
     85
     86  StatType *psfstats  = &results->psfstats;
     87  StatType *apstats   = &results->apstats;
    8188  StatType *kronstats = &results->kronstats;
    8289
    … …  
    125132    off_t meas = measureOffset;
    126133
    127     N = 0;
     134    int Nap = 0;
     135    int Npsf = 0;
     136    int Nkron = 0;
    128137    for (k = 0; k < averageT[0].Nmeasure; k++, meas++) {
    129138
    … …  
    164173      int myUbercalDist = getUbercalDist(meas, cat);
    165174      minUbercalDist = MIN(minUbercalDist, myUbercalDist);
     175
     176      int isUbercal = (measureT[k].dbFlags & ID_MEAS_PHOTOM_UBERCAL);
     177      float dMpsf = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR);
     178
     179      if (isUbercal) haveUbercal = TRUE; // haveUbercal is set per secfilt, isUbercal is per measure
    166180
    167181      if (isSetMrelFinal) {
    … …  
    231245          if (measureT[k].M < MaxMagForceSynth) {
    232246            forceSynth = TRUE;
    233             forceSynthEntry = N;
     247            forceSynthEntry = Npsf;
    234248          } else {
    235249            if (pass < 4) {
    … …  
    240254        }
    241255
    242         // Map (aplist) and Mkron (kronlist,dkronlist) are used to calculate mean mags per filter
     256        // Map (Maplist) and Mkron (Mkronlist,dkronlist) are used to calculate mean mags per filter
    243257        float Map = PhotAper (&measure[k]);
    244         aplist[N] = Map - Mcal - Mmos - Mgrid;
     258        if (!isnan(Map)) {
     259            Maplist[Nap] = Map - Mcal - Mmos - Mgrid;
     260            daplist[Nap] = dMpsf;
     261            waplist[Nap] = isUbercal ? 10.0 : 1.0;
     262            Nap ++;
     263        }
    245264
    246265        float Mkron = PhotKron (&measure[k]);
    247         kronlist[N] = Mkron - Mcal - Mmos - Mgrid;
    248         dkronlist[N] = measure[k].dMkron;
    249       }
     266        if (!isnan(Mkron)) {
     267          Mkronlist[Nkron] = Mkron - Mcal - Mmos - Mgrid;
     268          dkronlist[Nkron] = measure[k].dMkron;
     269          wkronlist[Nkron] = isUbercal ? 10.0 : 1.0;
     270          Nkron ++;
     271        }
     272      } // if (isSetMrelFinal)
    250273
    251274      // dlist gives the error per measurement, wlist gives the weight
    … …  
    255278      // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
    256279      // 4) some reference photcode of some kind can be specified as fixed and have a high weight
    257       dlist[N] = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR);
    258       wlist[N] = 1.0;
    259       list[N]  = Msys - Mcal - Mmos - Mgrid;
    260 
    261       // up-weight the ubercal values (or convergence can take a long time...)
    262       if (measureT[k].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
    263         haveUbercal = TRUE;
    264         wlist[N] = 10.0;
    265       }
     280      Mpsflist[Npsf] = Msys - Mcal - Mmos - Mgrid;
     281      dpsflist[Npsf] = dMpsf;
     282      wpsflist[Npsf] = isUbercal ? 10.0 : 1.0;
     283
    266284
    267285      // tie down reference photometry if the -refcode (code) option is selected
    268286      // eg, -refcode g_SDSS
    269287      // this probably makes no sense in the context of multifilter analysis
     288      // XXX probably need to use the photocde table to assign reference mag weights.
    270289      if (refPhotcode) {
    271290        if (code->code == refPhotcode->code) {
    272           wlist[N] = 100.0;
    273         }
    274       }
    275       N++;
     291          wpsflist[Npsf] = 100.0;
     292        }
     293      }
     294      Npsf ++;
    276295    }
    277296
    … …  
    288307
    289308    // when performing the grid analysis, STAR_TOOFEW should be set to 1;
    290     if (N < Nminmeas) { /* too few measurements */
    291       // fprintf (f, "%10.6f %10.6f %d %d %d\n", averageT[0].R, averageT[0].D, measureT[0].imageID, N, STAR_TOOFEW);
     309    if (Npsf < Nminmeas) { /* too few measurements */
     310      // fprintf (f, "%10.6f %10.6f %d %d %d\n", averageT[0].R, averageT[0].D, measureT[0].imageID, Npsf, STAR_TOOFEW);
    292311      secfilt[Nsec].flags |= ID_STAR_FEW;
    293312      if (Ncode == 0) {
    … …  
    303322    if (forceSynth) {
    304323      // use the single SYNTH value instead of the other mags here
    305       myAssert ((forceSynthEntry < N) && (forceSynthEntry >= 0), "programming error");
    306       list[0]  = list[forceSynthEntry];
    307       dlist[0] = dlist[forceSynthEntry];
    308       wlist[0] = wlist[forceSynthEntry];
    309       N = 1;
    310     }
    311     liststats (list, dlist, wlist, N, stats);
    312 
    313     secfilt[Nsec].M  = stats->mean;
    314     secfilt[Nsec].dM = stats->error;
    315     secfilt[Nsec].Xm = (stats->Nmeas > 1) ? 100.0*log10(stats->chisq + 1e-4) : NAN_S_SHORT;
     324      myAssert ((forceSynthEntry < Npsf) && (forceSynthEntry >= 0), "programming error");
     325      Mpsflist[0] = Mpsflist[forceSynthEntry];
     326      dpsflist[0] = dpsflist[forceSynthEntry];
     327      wpsflist[0] = wpsflist[forceSynthEntry];
     328      Npsf = 1;
     329    }
     330    liststats (Mpsflist, dpsflist, wpsflist, Npsf, psfstats);
     331
     332    secfilt[Nsec].M  = psfstats->mean;
     333    secfilt[Nsec].dM = psfstats->error;
     334    secfilt[Nsec].Xm = (psfstats->Nmeas > 1) ? 100.0*log10(psfstats->chisq + 1e-4) : NAN_S_SHORT;
    316335
    317336    // when running -averages, we have no information about the images, so we cannot set this
    … …  
    323342      found[Nsec] = TRUE;
    324343
    325       secfilt[Nsec].Mstdev = 1000.0*stats->sigma; // Mstdev is in millimags (not enough space for more precision)
     344      secfilt[Nsec].Mstdev = 1000.0*psfstats->sigma; // Mstdev is in millimags (not enough space for more precision)
    326345      // secfilt[Nsec].Ncode = Ncode;
    327       secfilt[Nsec].Nused = stats->Nmeas;
    328 
    329       secfilt[Nsec].M_80 = 1000 * stats->Upper80;
    330       secfilt[Nsec].M_20 = 1000 * stats->Lower20;
     346      secfilt[Nsec].Nused = psfstats->Nmeas;
     347
     348      secfilt[Nsec].M_80 = 1000 * psfstats->Upper80;
     349      secfilt[Nsec].M_20 = 1000 * psfstats->Lower20;
    331350
    332351      // NOTE : use the modified weight for apmags as well as psf mags
    333       liststats (aplist, dlist, wlist, N, apstats);
    334       secfilt[Nsec].Map  = apstats->mean;
    335 
    336       liststats (kronlist, dkronlist, wlist, N, kronstats);
    337       secfilt[Nsec].Mkron  = kronstats->mean;
    338       secfilt[Nsec].dMkron = kronstats->error;
     352      liststats (Maplist, daplist, waplist, Nap, apstats);
     353      secfilt[Nsec].Map  = Nap > 0 ? apstats->mean : NAN;
     354
     355      liststats (Mkronlist, dkronlist, wkronlist, Nkron, kronstats);
     356      secfilt[Nsec].Mkron  = Nkron > 0 ? kronstats->mean  : NAN;
     357      secfilt[Nsec].dMkron = Nkron > 0 ? kronstats->error : NAN;
    339358
    340359      // NOTE: for 2MASS measurements, Next should be 1, as should N
    341       if ((Next > 0) && (Next > 0.5*N)) {
     360      if ((Next > 0) && (Next > 0.5*Npsf)) {
    342361        secfilt[Nsec].flags |= ID_SECF_OBJ_EXT;
    343362      }
Note: See TracChangeset for help on using the changeset viewer.