IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36917


Ignore:
Timestamp:
Jun 19, 2014, 6:10:10 PM (12 years ago)
Author:
eugene
Message:

update to make relphot understand forced warps

Location:
branches/eam_branches/ipp-20140610/Ohana/src
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20140610/Ohana/src/libautocode/def/secfilt-ps1-v5.d

    r36899 r36917  
    22EXTNAME      DVO_SECFILT_PS1_V5
    33TYPE         BINTABLE
    4 SIZE         152
     4SIZE         160
    55DESCRIPTION  DVO SecFilt : Secondary Filter Data
    66
    … …  
    1010FIELD  Map,           MAG_AP,            float,    ave aperture mag in this band,         mags
    1111FIELD  dMap,          MAG_AP_ERR,        float,    ave aperture mag in this band,         mags
     12FIELD  sMap,          MAG_AP_STDEV,    float,      standard deviation of ap mags,         mags
    1213FIELD  Mkron,         MAG_KRON,          float,    ave kron mag in this band,             mags
    1314FIELD  dMkron,        MAG_KRON_ERR,      float,    formal error on average kron mag,      mags
     15FIELD  sMkron,        MAG_KRON_STDEV,  float,      standard deviation of kron mags,       mags
    1416                                         
    1517# these statistics are PSF-specific     

  • branches/eam_branches/ipp-20140610/Ohana/src/libautocode/def/secfilt.d

    r36899 r36917  
    22EXTNAME      DVO_SECFILT
    33TYPE         BINTABLE
    4 SIZE         152
     4SIZE         160
    55DESCRIPTION  DVO SecFilt : Secondary Filter Data
    66
    … …  
    1010FIELD  Map,           MAG_AP,          float,      average mag in this band,              mags
    1111FIELD  dMap,          MAG_AP_ERR,      float,      formal error on average mag,           mags
     12FIELD  sMap,          MAG_AP_STDEV,    float,      standard deviation of ap mags,         mags
    1213FIELD  Mkron,         MAG_KRON,        float,      ave kron mag in this band,             mags
    1314FIELD  dMkron,        MAG_KRON_ERR,    float,      formal error on average kron mag,      mags
     15FIELD  sMkron,        MAG_KRON_STDEV,  float,      standard deviation of kron mags,       mags
    1416
    1517# these statistics are PSF-specific     

  • branches/eam_branches/ipp-20140610/Ohana/src/libdvo/include/dvo.h

    r36912 r36917  
    599599  /* pointers for data manipulation */
    600600  off_t *found;
     601  off_t *foundWarp;
    601602  off_t *image;
    602603  off_t *mosaic;

  • branches/eam_branches/ipp-20140610/Ohana/src/libdvo/src/dvo_catalog.c

    r36898 r36917  
    166166  secfilt->Map         = NAN;
    167167  secfilt->dMap        = NAN;
     168  secfilt->sMap        = NAN;
    168169  secfilt->Mkron       = NAN;
    169170  secfilt->dMkron      = NAN;
     171  secfilt->sMkron      = NAN;
    170172
    171173  secfilt->Mstdev      = NAN;

  • branches/eam_branches/ipp-20140610/Ohana/src/libdvo/src/dvo_convert_PS1_V5.c

    r36900 r36917  
    284284    out[i].dM            = in[i].dM;     
    285285    out[i].Map           = in[i].Map;     
    286     out[i].dMap          = in[i].dM;     
     286    out[i].dMap          = in[i].dMap;     
     287    out[i].sMap          = in[i].sMap;     
    287288    out[i].Mkron         = in[i].Mkron;     
    288289    out[i].dMkron        = in[i].dMkron;     
    … …  
    329330    out[i].sFapWrp       = in[i].sFapWrp;
    330331
    331     out[i].NusedWrp      = in[i].Nused;
    332     out[i].NusedKronWrp  = in[i].NusedKron;
    333     out[i].NusedApWrp    = in[i].NusedAp;
     332    out[i].NusedWrp      = in[i].NusedWrp;
     333    out[i].NusedKronWrp  = in[i].NusedKronWrp;
     334    out[i].NusedApWrp    = in[i].NusedApWrp;
    334335
    335336    out[i].ubercalDist   = in[i].ubercalDist;     
    … …  
    350351    out[i].dM            = in[i].dM;     
    351352    out[i].Map           = in[i].Map;     
    352     out[i].dMap          = in[i].dM;     
     353    out[i].dMap          = in[i].dMap;     
     354    out[i].sMap          = in[i].sMap;     
    353355    out[i].Mkron         = in[i].Mkron;     
    354356    out[i].dMkron        = in[i].dMkron;     
     357    out[i].sMkron        = in[i].sMkron;     
    355358
    356359    out[i].Mstdev        = in[i].Mstdev;     
    … …  
    395398    out[i].sFapWrp       = in[i].sFapWrp;
    396399
    397     out[i].NusedWrp      = in[i].Nused;
    398     out[i].NusedKronWrp  = in[i].NusedKron;
    399     out[i].NusedApWrp    = in[i].NusedAp;
     400    out[i].NusedWrp      = in[i].NusedWrp;
     401    out[i].NusedKronWrp  = in[i].NusedKronWrp;
     402    out[i].NusedApWrp    = in[i].NusedApWrp;
    400403
    401404    out[i].ubercalDist   = in[i].ubercalDist;     

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

    r36910 r36917  
    153153  int i;
    154154
     155  fprintf (stderr, "this version of setMrel is invalid for now\n");
     156  exit (1);
     157       
    155158  int Nsecfilt = GetPhotcodeNsecfilt ();
    156159

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

    r36910 r36917  
    11# include "relphot.h"
     2# define UBERCAL_WEIGHT 100.0
    23
    34# define SKIP_THIS_MEAS(REASON) {                               \
    … …  
    1314    continue; }
    1415
     16static float MagToFlux (float Mag) {
     17  float Flux = pow(10.0, -0.4*(Mag));
     18  return (Flux);
     19}
     20
    1521int print_measure_set_alt (Average *average, SecFilt *secfilt, Measure *measure) {
    1622
    … …  
    2935  return (TRUE);
    3036}
     37
     38int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found);
    3139
    3240int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
    … …  
    5967      setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]);
    6068    }
     69
     70    // only measure force-warp mean values if issetMrelFinal (make it optional?)
     71    if (isSetMrelFinal) {
     72      setMrelAverageForcedWarp (m, Nc, pass, results, &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], &catalog[Nc].measure[m], &catalog[Nc].foundWarp[Nsecfilt*j]);
     73      setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]);
     74    }
    6175  }
    6276  if (primaryCell) free (primaryCell);
    … …  
    6478}
    6579
     80// set mean of chip measurements (selected by photcode range for now):
    6681int setMrelAverageExposure (off_t measureOffset, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found) {
    6782
    … …  
    199214        }
    200215
    201         // gpc1 stack data
    202         if ((measure[k].photcode >= 11000) && (measure[k].photcode <= 11400)) {
     216        // SKIP gpc1 stack data
     217        if ((measure[k].photcode >= 11000) && (measure[k].photcode <= 11500)) {
     218          continue;
     219        }
     220
     221        // SKIP gpc1 forced-warp data
     222        if ((measure[k].photcode >= 12000) && (measure[k].photcode <= 12500)) {
    203223          continue;
    204224        }
    … …  
    258278        if (!isnan(Map)) {
    259279            Maplist[Nap] = Map - Mcal - Mmos - Mgrid;
    260             daplist[Nap] = dMpsf;
    261             waplist[Nap] = isUbercal ? 10.0 : 1.0;
     280            daplist[Nap] = dMpsf;// XXX check on this...
     281            waplist[Nap] = isUbercal ? UBERCAL_WEIGHT : 1.0;
    262282            Nap ++;
    263283        }
    … …  
    267287          Mkronlist[Nkron] = Mkron - Mcal - Mmos - Mgrid;
    268288          dkronlist[Nkron] = measure[k].dMkron;
    269           wkronlist[Nkron] = isUbercal ? 10.0 : 1.0;
     289          wkronlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0;
    270290          Nkron ++;
    271291        }
    … …  
    280300      Mpsflist[Npsf] = Msys - Mcal - Mmos - Mgrid;
    281301      dpsflist[Npsf] = dMpsf;
    282       wpsflist[Npsf] = isUbercal ? 10.0 : 1.0;
     302      wpsflist[Npsf] = isUbercal ? UBERCAL_WEIGHT : 1.0;
    283303
    284304      // NOTE:
    … …  
    293313      if (refPhotcode) {
    294314        if (code->code == refPhotcode->code) {
    295           wpsflist[Npsf] = 100.0;
     315          wpsflist[Npsf] = UBERCAL_WEIGHT;
    296316        }
    297317      }
    … …  
    304324    if (isSetMrelFinal && (pass == 0)) {
    305325      if ((thisCode < 6) || (thisCode == 9)) {
    306         secfilt[Nsec].Ncode = NexpPS1;
     326        secfilt[Nsec].Ncode = NexpPS1; 
    307327      } else {
    308328        secfilt[Nsec].Ncode = Ncode; // 2MASS data if it exists
    … …  
    346366      found[Nsec] = TRUE;
    347367
    348       secfilt[Nsec].Mstdev = 1000.0*psfstats->sigma; // Mstdev is in millimags (not enough space for more precision)
     368      secfilt[Nsec].Mstdev = psfstats->sigma; // Mstdev is in millimags (not enough space for more precision)
    349369      // secfilt[Nsec].Ncode = Ncode;
    350370      secfilt[Nsec].Nused = psfstats->Nmeas;
    351371
    352       secfilt[Nsec].Mmax = 1000 * psfstats->max;
    353       secfilt[Nsec].Mmin = 1000 * psfstats->min;
     372      secfilt[Nsec].Mmax = psfstats->max;
     373      secfilt[Nsec].Mmin = psfstats->min;
    354374
    355375      // NOTE : use the modified weight for apmags as well as psf mags
    356376      liststats (Maplist, daplist, waplist, Nap, apstats);
    357377      secfilt[Nsec].Map  = Nap > 0 ? apstats->mean : NAN;
     378      secfilt[Nsec].dMap  = Nap > 0 ? apstats->error : NAN;
     379      secfilt[Nsec].sMap  = Nap > 0 ? apstats->sigma : NAN;
     380      secfilt[Nsec].NusedAp  = Nap;
    358381
    359382      liststats (Mkronlist, dkronlist, wkronlist, Nkron, kronstats);
    360383      secfilt[Nsec].Mkron  = Nkron > 0 ? kronstats->mean  : NAN;
    361384      secfilt[Nsec].dMkron = Nkron > 0 ? kronstats->error : NAN;
     385      secfilt[Nsec].sMkron = Nkron > 0 ? kronstats->sigma : NAN;
     386      secfilt[Nsec].NusedKron  = Nkron;
    362387
    363388      // NOTE: for 2MASS measurements, Next should be 1, as should N
    … …  
    585610    secfilt[Nsec].FkronStk  = zpFactor * measure[k].FluxKron;
    586611    secfilt[Nsec].dFkronStk = zpFactor * measure[k].dFluxKron;
     612    secfilt[Nsec].FapStk    = zpFactor * measure[k].FluxAp;
     613    secfilt[Nsec].dFapStk   = zpFactor * measure[k].dFluxAp;
     614
     615    secfilt[Nsec].MpsfStk   = (measure[k].FluxPSF  > 0.0) ? zp -2.5*log10(measure[k].FluxPSF) : NAN;
     616    secfilt[Nsec].MkronStk  = (measure[k].FluxKron > 0.0) ? zp -2.5*log10(measure[k].FluxKron) : NAN;
     617    secfilt[Nsec].MapStk    = (measure[k].FluxAp   > 0.0) ? zp -2.5*log10(measure[k].FluxAp) : NAN;
    587618
    588619    secfilt[Nsec].stackDetectID = ID;
    … …  
    607638}
    608639
     640# undef SKIP_THIS_MEAS
     641# define SKIP_THIS_MEAS(REASON) {                               \
     642    measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM;  \
     643    results->REASON ++;                                         \
     644    continue; }
     645
     646// set mean of forced-warp measurements (selected by photcode range for now):
     647// somewhat simplified relative to chip-photometry:
     648// * no grid, no mosaic, no 2MASS, no SYNTH, no Ubercal, no flatcorr
     649int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found) {
     650
     651  off_t k;
     652  float Fsys = 0, Mcal= 0;
     653
     654  // we are measuring means for 3 types of FLUXes: psf, ap, kron.  I am using the psf mag
     655  // error for the ap mags, but krons have their own errors.  it is an open question if I
     656  // should be doing weighted or unweighted fits (this is a user option)
     657  double *Fpsflist  = results->Mpsflist;
     658  double *dpsflist  = results->dpsflist;
     659  double *wpsflist  = results->wpsflist;
     660
     661  double *Faplist   = results->Maplist;
     662  double *daplist   = results->daplist;
     663  double *waplist   = results->waplist;
     664
     665  double *Fkronlist = results->Mkronlist;
     666  double *dkronlist = results->dkronlist;
     667  double *wkronlist = results->wkronlist;
     668
     669  StatType *psfstats  = &results->psfstats;
     670  StatType *apstats   = &results->apstats;
     671  StatType *kronstats = &results->kronstats;
     672
     673  // option for a test print
     674  if (FALSE && (average[0].objID == 0x7146) && (average[0].catID == 0x49d8)) {
     675    fprintf (stderr, "test obj\n");
     676    print_measure_set_alt (average, secfilt, measure);
     677  }
     678
     679  int Ns;
     680  for (Ns = 0; Ns < Nphotcodes; Ns++) {
     681
     682    int thisCode = photcodes[Ns][0].code;
     683    int Nsec = GetPhotcodeNsec(thisCode);
     684
     685    /* calculate the average mag in this SEC photcode for a single star */
     686
     687    /* star/photcodes already calibrated */
     688    if (found[Nsec]) continue; 
     689     
     690    off_t meas = measureOffset;
     691
     692    int Nap = 0;
     693    int Npsf = 0;
     694    int Nkron = 0;
     695    for (k = 0; k < average[0].Nmeasure; k++, meas++) {
     696
     697      // skip measurements which are not ForcedWarp values
     698      if ((measure[k].photcode < 12000) || (measure[k].photcode > 12500)) continue;
     699
     700      // skip measurements that do not match the current photcode
     701      PhotCode *code = GetPhotcodebyCode (measure[k].photcode);
     702      if (!code) continue;
     703      if (code->equiv != thisCode) { continue; }
     704
     705      if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad);
     706
     707      if (getImageEntry (meas, cat) < 0) {
     708        // measurements without an image are either external reference photometry or
     709        // data for which the associated image has not been loaded (probably because of
     710        // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
     711        Mcal = measure[k].Mcal; // check that this is zero for loaded REF value
     712      } else {
     713        // use getMcal not getMcal_alt?
     714        Mcal  = getMcal_alt (meas, cat, NULL, measure[k].Xccd, measure[k].Yccd);
     715        // Mcal  = getMcal (meas, cat);
     716        if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
     717      }
     718      float Fcal = MagToFlux(Mcal);
     719
     720      // in the calculations below,
     721      // ...list gives the error per measurement, wlist gives the weight
     722      // we can modify the error and weight in a few ways:
     723      // 1) MIN_ERROR guarantees a floor
     724      // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
     725
     726      // skip some absurd values NAN
     727      // NOTE : I am using PhotCat not PhotSys for now since GPC1 chip-to-chip color terms
     728      // are small (and not measured)
     729      float Fpsf = PhotFluxCat (&measure[k], MAG_CLASS_PSF);
     730      if (isnan(Fsys)) SKIP_THIS_MEAS(Nsys);
     731      // if (Msys <  0.0) SKIP_THIS_MEAS(Nsys);
     732      // if (Msys > 30.0) SKIP_THIS_MEAS(Nsys);
     733
     734      float dFpsf = PhotFluxCatErr (&measure[k], MAG_CLASS_PSF);
     735      dFpsf = MAX (dFpsf, MIN_ERROR*Fpsf); // MIN_ERROR is a fractional error
     736      Fpsflist[Npsf] = Fpsf * Fcal;
     737      dpsflist[Npsf] = dFpsf * Fcal;
     738      wpsflist[Npsf] = 1.0;
     739      Npsf ++;
     740
     741      float Fap = PhotFluxCat (&measure[k], MAG_CLASS_APER);
     742      float dFap = PhotFluxCatErr (&measure[k], MAG_CLASS_APER);
     743      if (!isnan(Fap)) {
     744        Faplist[Nap] = Fap * Fcal;
     745        daplist[Nap] = dFap * Fcal;
     746        waplist[Nap] = 1.0; // drop weight lists?
     747        Nap ++;
     748      }
     749
     750      float Fkron = PhotFluxCat (&measure[k], MAG_CLASS_KRON);
     751      float dFkron = PhotFluxCatErr (&measure[k], MAG_CLASS_KRON);
     752      if (!isnan(Fkron)) {
     753        Fkronlist[Nkron] = Fkron * Fcal;
     754        dkronlist[Nkron] = dFkron * Fcal;
     755        wkronlist[Nkron] = 1.0;
     756        Nkron ++;
     757      }
     758    }
     759    if (Npsf < 1) continue;
     760
     761    found[Nsec] = TRUE;
     762
     763    liststats (Fpsflist, dpsflist, wpsflist, Npsf, psfstats);
     764
     765    secfilt[Nsec].FpsfWrp  = psfstats->mean;
     766    secfilt[Nsec].dFpsfWrp = psfstats->error;
     767    secfilt[Nsec].sFpsfWrp = psfstats->sigma; // Mstdev is in millimags (not enough space for more precision)
     768    secfilt[Nsec].NusedWrp = psfstats->Nmeas;
     769    secfilt[Nsec].MpsfWrp  = isnan(secfilt[Nsec].FpsfWrp) ? NAN : -2.5*log(secfilt[Nsec].FpsfWrp); // XXX need to add zero-point
     770
     771    // NOTE : use the modified weight for apmags as well as psf mags
     772    liststats (Faplist, daplist, waplist, Nap, apstats);
     773    secfilt[Nsec].FapWrp     = Nap > 0 ? apstats->mean : NAN;
     774    secfilt[Nsec].dFapWrp    = Nap > 0 ? apstats->error : NAN;
     775    secfilt[Nsec].sFapWrp    = Nap > 0 ? apstats->sigma : NAN;
     776    secfilt[Nsec].NusedApWrp = Nap;
     777    secfilt[Nsec].MapWrp  = isnan(secfilt[Nsec].FapWrp) ? NAN : -2.5*log(secfilt[Nsec].FapWrp); // XXX need to add zero-point
     778
     779    liststats (Fkronlist, dkronlist, wkronlist, Nkron, kronstats);
     780    secfilt[Nsec].FkronWrp     = Nkron > 0 ? kronstats->mean  : NAN;
     781    secfilt[Nsec].dFkronWrp    = Nkron > 0 ? kronstats->error : NAN;
     782    secfilt[Nsec].sFkronWrp    = Nkron > 0 ? kronstats->sigma : NAN;
     783    secfilt[Nsec].NusedKronWrp = Nkron;
     784    secfilt[Nsec].MkronWrp  = isnan(secfilt[Nsec].FkronWrp) ? NAN : -2.5*log(secfilt[Nsec].FkronWrp); // XXX need to add zero-point
     785  }
     786  return (TRUE);
     787}
     788
    609789int setGlobalObjStats (Average *average, Measure *measure) {
    610790

  • branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/setMrelFinal.c

    r36910 r36917  
    8686  /* set catalog[0].found[i] = FALSE */
    8787  ALLOCATE (catalog[0].found, off_t, MAX (1, Nsecfilt*catalog[0].Naverage));
     88  ALLOCATE (catalog[0].foundWarp, off_t, MAX (1, Nsecfilt*catalog[0].Naverage));
    8889  for (i = 0; i < Nsecfilt*catalog[0].Naverage; i++) {
    8990    catalog[0].found[i] = FALSE;
     91    catalog[0].foundWarp[i] = FALSE;
    9092  }
    9193
Note: See TracChangeset for help on using the changeset viewer.