IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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

merge from trunk

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

Legend:

Unmodified
Added
Removed

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

    r36520 r36535  
    9595  // others are only used on the final output steps.
    9696
    97   double *list;               // list of measure.mag values for a given star
    98   double *dlist;              // mag errors for a star
    99   double *wlist;              // weights to use for mean mags
    100   double *aplist;             // ap mags for a star
    101   double *kronlist;           // kron mags for a star
     97  double *Mpsflist;           // list of measure.mag values for a given star
     98  double *dpsflist;           // mag errors for a star
     99  double *wpsflist;           // weights to use for mean mags
     100
     101  double *Maplist;            // ap mags for a star
     102  double *daplist;            // ap mags for a star
     103  double *waplist;            // ap mags for a star
     104
     105  double *Mkronlist;          // kron mags for a star
    102106  double *dkronlist;          // kron mag errors
     107  double *wkronlist;          // kron mag errors
    103108
    104109  double *psfqf_list;         // psfqf for all filters
    … …  
    106111  double *stargal_list;       // stargal for all filters
    107112
    108   StatType stats;
     113  StatType psfstats;
    109114  StatType apstats;
    110115  StatType kronstats;

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

    r36519 r36535  
    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-20140206/Ohana/src/relphot/src/setMrelCatalog.c

    r36503 r36535  
    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;
     280      Mpsflist[Npsf] = Msys - Mcal - Mmos - Mgrid;
     281      dpsflist[Npsf] = dMpsf;
     282      wpsflist[Npsf] = isUbercal ? 10.0 : 1.0;
    260283
    261284      // NOTE:
    … …  
    264287      // Mmos and Mgrid are offsets for mosaic and grid
    265288
    266       // up-weight the ubercal values (or convergence can take a long time...)
    267       if (measureT[k].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
    268         haveUbercal = TRUE;
    269         wlist[N] = 10.0;
    270       }
    271 
    272289      // tie down reference photometry if the -refcode (code) option is selected
    273290      // eg, -refcode g_SDSS
    274291      // this probably makes no sense in the context of multifilter analysis
     292      // XXX probably need to use the photocde table to assign reference mag weights.
    275293      if (refPhotcode) {
    276294        if (code->code == refPhotcode->code) {
    277           wlist[N] = 100.0;
    278         }
    279       }
    280       N++;
     295          wpsflist[Npsf] = 100.0;
     296        }
     297      }
     298      Npsf ++;
    281299    }
    282300
    … …  
    293311
    294312    // when performing the grid analysis, STAR_TOOFEW should be set to 1;
    295     if (N < Nminmeas) { /* too few measurements */
    296       // fprintf (f, "%10.6f %10.6f %d %d %d\n", averageT[0].R, averageT[0].D, measureT[0].imageID, N, STAR_TOOFEW);
     313    if (Npsf < Nminmeas) { /* too few measurements */
     314      // fprintf (f, "%10.6f %10.6f %d %d %d\n", averageT[0].R, averageT[0].D, measureT[0].imageID, Npsf, STAR_TOOFEW);
    297315      secfilt[Nsec].flags |= ID_STAR_FEW;
    298316      if (Ncode == 0) {
    … …  
    308326    if (forceSynth) {
    309327      // use the single SYNTH value instead of the other mags here
    310       myAssert ((forceSynthEntry < N) && (forceSynthEntry >= 0), "programming error");
    311       list[0]  = list[forceSynthEntry];
    312       dlist[0] = dlist[forceSynthEntry];
    313       wlist[0] = wlist[forceSynthEntry];
    314       N = 1;
    315     }
    316     liststats (list, dlist, wlist, N, stats);
    317 
    318     secfilt[Nsec].M  = stats->mean;
    319     secfilt[Nsec].dM = stats->error;
    320     secfilt[Nsec].Xm = (stats->Nmeas > 1) ? 100.0*log10(stats->chisq + 1e-4) : NAN_S_SHORT;
     328      myAssert ((forceSynthEntry < Npsf) && (forceSynthEntry >= 0), "programming error");
     329      Mpsflist[0] = Mpsflist[forceSynthEntry];
     330      dpsflist[0] = dpsflist[forceSynthEntry];
     331      wpsflist[0] = wpsflist[forceSynthEntry];
     332      Npsf = 1;
     333    }
     334    liststats (Mpsflist, dpsflist, wpsflist, Npsf, psfstats);
     335
     336    secfilt[Nsec].M  = psfstats->mean;
     337    secfilt[Nsec].dM = psfstats->error;
     338    secfilt[Nsec].Xm = (psfstats->Nmeas > 1) ? 100.0*log10(psfstats->chisq + 1e-4) : NAN_S_SHORT;
    321339
    322340    // when running -averages, we have no information about the images, so we cannot set this
    … …  
    328346      found[Nsec] = TRUE;
    329347
    330       secfilt[Nsec].Mstdev = 1000.0*stats->sigma; // Mstdev is in millimags (not enough space for more precision)
     348      secfilt[Nsec].Mstdev = 1000.0*psfstats->sigma; // Mstdev is in millimags (not enough space for more precision)
    331349      // secfilt[Nsec].Ncode = Ncode;
    332       secfilt[Nsec].Nused = stats->Nmeas;
    333 
    334       secfilt[Nsec].M_80 = 1000 * stats->Upper80;
    335       secfilt[Nsec].M_20 = 1000 * stats->Lower20;
     350      secfilt[Nsec].Nused = psfstats->Nmeas;
     351
     352      secfilt[Nsec].M_80 = 1000 * psfstats->Upper80;
     353      secfilt[Nsec].M_20 = 1000 * psfstats->Lower20;
    336354
    337355      // NOTE : use the modified weight for apmags as well as psf mags
    338       liststats (aplist, dlist, wlist, N, apstats);
    339       secfilt[Nsec].Map  = apstats->mean;
    340 
    341       liststats (kronlist, dkronlist, wlist, N, kronstats);
    342       secfilt[Nsec].Mkron  = kronstats->mean;
    343       secfilt[Nsec].dMkron = kronstats->error;
     356      liststats (Maplist, daplist, waplist, Nap, apstats);
     357      secfilt[Nsec].Map  = Nap > 0 ? apstats->mean : NAN;
     358
     359      liststats (Mkronlist, dkronlist, wkronlist, Nkron, kronstats);
     360      secfilt[Nsec].Mkron  = Nkron > 0 ? kronstats->mean  : NAN;
     361      secfilt[Nsec].dMkron = Nkron > 0 ? kronstats->error : NAN;
    344362
    345363      // NOTE: for 2MASS measurements, Next should be 1, as should N
    346       if ((Next > 0) && (Next > 0.5*N)) {
     364      if ((Next > 0) && (Next > 0.5*Npsf)) {
    347365        secfilt[Nsec].flags |= ID_SECF_OBJ_EXT;
    348366      }
Note: See TracChangeset for help on using the changeset viewer.