Changeset 38289

Timestamp:

May 18, 2015, 8:26:24 PM (11 years ago)

Author:

eugene

Message:

changing setMrelCatalog logic

File:

• branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelCatalog.c (modified) (3 diffs)

branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelCatalog.c

@@ -86 +86 @@
   // error for the ap mags, but krons have their own errors.  it is an open question if I
   // should be doing weighted or unweighted fits (this is a user option)
-  double *Mpsflist  = results->Mpsflist;
-  double *dpsflist  = results->dpsflist;
-  double *wpsflist  = results->wpsflist;
-
-  double *Maplist   = results->Maplist;
-  double *daplist   = results->daplist;
-  double *waplist   = results->waplist;
-
-  double *Mkronlist = results->Mkronlist;
-  double *dkronlist = results->dkronlist;
-  double *wkronlist = results->wkronlist;
+
+  // these lists are pre-allocated to the max length
+  double **Mpsflist  = results->Mpsflist;
+  double **dpsflist  = results->dpsflist;
+  double **wpsflist  = results->wpsflist;
+
+  double **Maplist   = results->Maplist;
+  double **daplist   = results->daplist;
+  double **waplist   = results->waplist;
+
+  double **Mkronlist = results->Mkronlist;
+  double **dkronlist = results->dkronlist;
+  double **wkronlist = results->wkronlist;
 
   StatType *psfstats  = &results->psfstats;
@@ -108 +110 @@
     fprintf (stderr, "test obj\n");
     print_measure_set_alt (average, secfilt, measure);
+  }
+
+  // first pass: check for PS1[n], Synth, Tycho, 2MASS
+
+  // first, determine the state of things for each of the main filters:
+  // in this bit, we are using hardwired photcode values:
+  // PS1_g, PS1_r, PS1_i, PS1_z, PS1_y, PS1_w == (1,2,3,4,5,6); 0 == none
+  // note that Nsec = PS_N - 1
+
+  int havePS1[7];
+  int haveSYN[7];
+  int needSYN[7];
+  int measSYN[7];
+
+  int haveTYCHO = FALSE;
+  int have2MASS = FALSE;
+
+  // was: 13.64, 13.76, 13.74, 12.94, 12.01
+  float MaxMagForceSynth[] = {0.0, 13.5, 13.5, 13.5, 13.0, 12.0, 13.5};
+  //                          NONE, g,    r,    i,    z,    y,    w
+
+  for (i = 0; i < 7; i++) {
+    havePS1[i] = FALSE; // do we have any PS1 measurement for this filter?
+    haveSYN[i] = FALSE; // do we have any synthetic measurement for this filter?
+    needSYN[i] = FALSE; // force the use of the synthetic measurement for this filter?
+    measSYN[i] = -1; // which entries carries the synthetic measurement for this filter?
+  }
+
+  for (k = 0; k < average[0].Nmeasure; k++) {
+
+    // skip measurements that do not have a valid photcode (which ones lack photcode?)
+    PhotCode *code = GetPhotcodebyCode (measure[k].photcode);
+    if (!code) continue;
+
+    // is this synth_grizy?
+    // is this PS1?
+
+    // are we a PS1 exposure photcode?
+    if (isGPC1chip(measure[k].photcode)) {
+      int Nfilter = whichGPC1filter (measure[k].photcode);
+      havePS1[Nfilter] = TRUE;
+    }
+    if (isGPC1synth(measure[k].photcode)) {
+      int Nfilter = measure[k].photcode - 3000;
+      haveSYN[Nfilter] = TRUE;
+      measSYN[Nfilter] = k;
+
+      // force the use of SYN even if we have PS1 mags?
+      if (measure[k].M < MaxMagForceSynth[Nfilter]) {
+        needSYN[Nfilter] = TRUE;
+      } 
+    }
+    if (isTycho(measure[k].photcode)) {
+      haveTYCHO = TRUE;
+    }
+    if (is2MASS(measure[k].photcode)) {
+      have2MASS = TRUE;
+    }
   }
 
@@ -450 +510 @@
 }
 
+// set mean of chip measurements (selected by photcode range for now):
+int setMrelAverageExposure_old (off_t measureOffset, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found) {
+
+  off_t k;
+  float Msys = 0, Mcal= 0, Mmos = 0, Mgrid = 0;
+
+  // we are measuring means for 3 types of mags: psf, ap, kron.  I am using the psf mag
+  // error for the ap mags, but krons have their own errors.  it is an open question if I
+  // should be doing weighted or unweighted fits (this is a user option)
+  double *Mpsflist  = results->Mpsflist;
+  double *dpsflist  = results->dpsflist;
+  double *wpsflist  = results->wpsflist;
+
+  double *Maplist   = results->Maplist;
+  double *daplist   = results->daplist;
+  double *waplist   = results->waplist;
+
+  double *Mkronlist = results->Mkronlist;
+  double *dkronlist = results->dkronlist;
+  double *wkronlist = results->wkronlist;
+
+  StatType *psfstats  = &results->psfstats;
+  StatType *apstats   = &results->apstats;
+  StatType *kronstats = &results->kronstats;
+
+  int isSetMrelFinal = (pass >= 0);
+
+  // option for a test print
+  if (FALSE && (average[0].objID == 0x7146) && (average[0].catID == 0x49d8)) {
+    fprintf (stderr, "test obj\n");
+    print_measure_set_alt (average, secfilt, measure);
+  }
+
+  int GoodPS1 = FALSE;
+  int Good2MASS = FALSE;
+  int Galaxy2MASS = FALSE;
+
+  int NextPS1 = 0;
+  int NpsfPS1 = 0;
+
+  int Ns;
+  for (Ns = 0; Ns < Nphotcodes; Ns++) {
+
+    int thisCode = photcodes[Ns][0].code;
+    int Nsec = GetPhotcodeNsec(thisCode);
+
+    /* calculate the average mag in this SEC photcode for a single star */
+
+    /* star/photcodes already calibrated */
+    if ( isSetMrelFinal && found[Nsec]) continue;  
+      
+    // skip bad stars
+    if (!isSetMrelFinal && (secfilt[Nsec].flags & STAR_BAD)) continue;
+
+    int NexpPS1 = 0;
+    int Ncode = 0;
+    int Next = 0;
+    int haveSynth = FALSE;
+    int havePS1   = FALSE;
+
+    int forceSynth = FALSE;
+    int forceSynthEntry = -1;
+
+    int haveUbercal = FALSE;
+
+    int minUbercalDist = 1000;
+    
+    off_t meas = measureOffset;
+
+    float psfQfMax = 0.0;
+    float psfQfPerfMax = 0.0;
+
+    int Nap = 0;
+    int Npsf = 0;
+    int Nkron = 0;
+    for (k = 0; k < averageT[0].Nmeasure; k++, meas++) {
+
+      // skip measurements that do not match the current photcode
+      PhotCode *code = GetPhotcodebyCode (measureT[k].photcode);
+      if (!code) continue;
+      if (code->equiv != thisCode) { continue; }
+      Ncode ++;
+
+      // are we a PS1 exposure photcode?
+      if (isGPC1chip(measureT[k].photcode)) NexpPS1 ++;
+
+      // SKIP gpc1 stack data
+      if (isGPC1stack(measureT[k].photcode)) continue;
+
+      // SKIP gpc1 forced-warp data
+      if (isGPC1warp(measureT[k].photcode)) continue;
+
+      if (isSetMrelFinal) {
+        if (measure[k].psfQF     > psfQfMax)     psfQfMax     = measure[k].psfQF;
+        if (measure[k].psfQFperf > psfQfPerfMax) psfQfPerfMax = measure[k].psfQFperf;
+      }
+
+      if (measureT[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); 
+
+      if (getImageEntry (meas, cat) < 0) {
+        // measurements without an image are either external reference photometry or
+        // data for which the associated image has not been loaded (probably because of
+        // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
+        Mmos = Mgrid = 0;
+        Mcal = measureT[k].Mcal; // check that this is zero for loaded REF value
+      } else {
+        Mcal  = getMcal_alt (meas, cat, flatcorr, measureT[k].Xccd, measureT[k].Yccd);
+        if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
+        Mmos  = getMmos  (meas, cat);
+        if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
+        Mgrid = getMgrid (meas, cat);
+        if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
+      }
+
+      // skip some absurd values NAN, < 0.0, > 30.0
+      Msys = PhotSysTiny (&measureT[k], &averageT[0], &secfilt[0], MAG_CLASS_PSF);
+      if (isnan(Msys)) SKIP_THIS_MEAS(Nsys);
+      if (Msys <  0.0) SKIP_THIS_MEAS(Nsys);
+      if (Msys > 30.0) SKIP_THIS_MEAS(Nsys);
+
+      int myUbercalDist = getUbercalDist(meas, cat);
+      minUbercalDist = MIN(minUbercalDist, myUbercalDist);
+
+      int isUbercal = (measureT[k].dbFlags & ID_MEAS_PHOTOM_UBERCAL);
+      float dMpsf = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR);
+
+      if (isUbercal) haveUbercal = TRUE; // haveUbercal is set per secfilt, isUbercal is per measure
+
+      if (isSetMrelFinal) {
+        // special options for PS1 data
+        if ((measure[k].photcode >= 10000) && (measure[k].photcode <= 10500)) {
+          // count the extended detections
+          if (!isnan(measure[k].Map)) {
+            float dMagAp = measure[k].M - measure[k].Map;
+            float SigmaAp = hypot(0.1, 2.5*measure[k].dM);
+            // XXX this is still quite ad hoc, but at least it:
+            // (a) converges to 0.1 mag offset at the bright end
+            // (b) converges to 0.5 mag offset at the faint end (dM = 0.2)
+            if (dMagAp > SigmaAp) {
+              Next ++;
+              NextPS1 ++;
+            } else {
+              NpsfPS1 ++;
+            }
+          }
+          havePS1 = TRUE;
+        }
+
+        // count extended detections for 2MASS (XXX NOTE hardwired photcodes 2011, 2012, 2013)
+        if ((measure[k].photcode >= 2011) && (measure[k].photcode <= 2013)) {
+          if (measure[k].photFlags & 0x00c00000) {
+            Next ++;
+            Galaxy2MASS = TRUE;
+          }
+          if (pass == 0) {
+            if (measure[k].photFlags & 0x00000007) {
+              Good2MASS = TRUE;
+            } else {
+              // detections without one of these bits should only be used in PASS_1
+              SKIP_THIS_MEAS(Nbad);
+            }
+          }
+        }
+
+        // Blindly accepth the SYNTH mags if we are above saturation, otherwise, 
+        // ignore SYNTH photcodes until PASS == 4 (where we also accept saturated stars)
+        if ((measure[k].photcode >= 3001) && (measure[k].photcode <= 3005)) {
+          // something of a hack: force object to use synth values if synth mags >>
+          // saturation (3pi instrumental mags < -15)
+          float MaxMagForceSynth = NAN;
+          switch (measure[k].photcode) {
+            case 3001:
+              MaxMagForceSynth = 13.64;
+              break;
+            case 3002:
+              MaxMagForceSynth = 13.76;
+              break;
+            case 3003:
+              MaxMagForceSynth = 13.74;
+              break;
+            case 3004:
+              MaxMagForceSynth = 12.94;
+              break;
+            case 3005:
+              MaxMagForceSynth = 12.01;
+              break;
+          }
+          if (measureT[k].M < MaxMagForceSynth) {
+            forceSynth = TRUE;
+            forceSynthEntry = Npsf;
+          } else {
+            if (pass < 4) {
+              SKIP_THIS_MEAS(Nbad);
+            }
+            haveSynth = TRUE;
+          }
+        }
+
+        // Map (Maplist) and Mkron (Mkronlist,dkronlist) are used to calculate mean mags per filter
+        float Map = PhotCat (&measure[k], MAG_CLASS_APER);
+        if (!isnan(Map)) {
+            Maplist[Nap] = Map - Mcal - Mmos - Mgrid;
+            daplist[Nap] = dMpsf;// XXX check on this...
+            waplist[Nap] = isUbercal ? UBERCAL_WEIGHT : 1.0;
+            Nap ++;
+        }
+
+        float Mkron = PhotCat (&measure[k], MAG_CLASS_KRON);
+        if (!isnan(Mkron)) {
+          Mkronlist[Nkron] = Mkron - Mcal - Mmos - Mgrid;
+          dkronlist[Nkron] = measure[k].dMkron;
+          wkronlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0;
+          Nkron ++;
+        }
+      } // if (isSetMrelFinal)
+
+      // dlist gives the error per measurement, wlist gives the weight
+      // we can modify the error and weight in a few ways:
+      // 1) MIN_ERROR guarantees a floor
+      // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
+      // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
+      // 4) some reference photcode of some kind can be specified as fixed and have a high weight
+      Mpsflist[Npsf] = Msys - Mcal - Mmos - Mgrid;
+      dpsflist[Npsf] = dMpsf;
+      wpsflist[Npsf] = isUbercal ? UBERCAL_WEIGHT : 1.0;
+
+      // NOTE: 
+      // Msys is measure[i].M + zp corrections
+      // Mcal is image[j].Mcal
+      // Mmos and Mgrid are offsets for mosaic and grid
+
+      // tie down reference photometry if the -refcode (code) option is selected
+      // eg, -refcode g_SDSS
+      // this probably makes no sense in the context of multifilter analysis
+      // XXX probably need to use the photocde table to assign reference mag weights.
+      if (refPhotcode) {
+        if (code->code == refPhotcode->code) {
+          wpsflist[Npsf] = UBERCAL_WEIGHT;
+        }
+      }
+      Npsf ++;
+    }
+
+    int Nminmeas = isSetMrelFinal ? 1 : STAR_TOOFEW + 1;
+
+    // XXX : ugh : another hard-wired photcode entry...
+    if (isSetMrelFinal && (pass == 0)) {
+      if ((thisCode < 6) || (thisCode == 9)) {
+        secfilt[Nsec].Ncode = NexpPS1; 
+      } else {
+        secfilt[Nsec].Ncode = Ncode; // 2MASS data if it exists
+      }
+    }
+
+    // when performing the grid analysis, STAR_TOOFEW should be set to 1;
+    if (Npsf < Nminmeas) { /* too few measurements */
+      // fprintf (f, "%10.6f %10.6f %d %d %d\n", averageT[0].R, averageT[0].D, measureT[0].imageID, Npsf, STAR_TOOFEW); 
+      secfilt[Nsec].flags |= ID_STAR_FEW;
+      if (Ncode == 0) {
+        results->Ncode ++;
+      } else {
+        results->Nfew ++;
+      }
+      continue;
+    } else {
+      secfilt[Nsec].flags &= ~ID_STAR_FEW;
+    }   
+
+    if (forceSynth) {
+      // use the single SYNTH value instead of the other mags here
+      myAssert ((forceSynthEntry < Npsf) && (forceSynthEntry >= 0), "programming error");
+      Mpsflist[0] = Mpsflist[forceSynthEntry];
+      dpsflist[0] = dpsflist[forceSynthEntry];
+      wpsflist[0] = wpsflist[forceSynthEntry];
+      Npsf = 1;
+    }
+    liststats (Mpsflist, dpsflist, wpsflist, Npsf, psfstats);
+
+    secfilt[Nsec].M      = psfstats->mean;
+    secfilt[Nsec].dM     = psfstats->error;
+    secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
+
+    // when running -averages, we have no information about the images, so we cannot set this
+    if (minUbercalDist > -1) {
+      secfilt[Nsec].ubercalDist = minUbercalDist;
+    }
+
+    if (isSetMrelFinal) {
+      found[Nsec] = TRUE;
+
+      secfilt[Nsec].Mstdev = psfstats->sigma; // Mstdev is in millimags (not enough space for more precision)
+      // secfilt[Nsec].Ncode = Ncode;
+      secfilt[Nsec].Nused = psfstats->Nmeas;
+
+      secfilt[Nsec].Mmax = psfstats->max;
+      secfilt[Nsec].Mmin = psfstats->min;
+
+      secfilt[Nsec].psfQfMax     = psfQfMax;
+      secfilt[Nsec].psfQfPerfMax = psfQfPerfMax;
+
+      // NOTE : use the modified weight for apmags as well as psf mags
+      liststats (Maplist, daplist, waplist, Nap, apstats);
+      secfilt[Nsec].Map  = Nap > 0 ? apstats->mean : NAN; 
+      secfilt[Nsec].dMap  = Nap > 0 ? apstats->error : NAN; 
+      secfilt[Nsec].sMap  = Nap > 0 ? apstats->sigma : NAN; 
+      secfilt[Nsec].NusedAp  = Nap;
+
+      liststats (Mkronlist, dkronlist, wkronlist, Nkron, kronstats);
+      secfilt[Nsec].Mkron  = Nkron > 0 ? kronstats->mean  : NAN; 
+      secfilt[Nsec].dMkron = Nkron > 0 ? kronstats->error : NAN; 
+      secfilt[Nsec].sMkron = Nkron > 0 ? kronstats->sigma : NAN; 
+      secfilt[Nsec].NusedKron  = Nkron;
+
+      // NOTE: for 2MASS measurements, Next should be 1, as should N
+      if ((Next > 0) && (Next > 0.5*Npsf)) {
+        secfilt[Nsec].flags |= ID_SECF_OBJ_EXT;
+      }
+
+      switch (pass) {
+        case 0:
+          secfilt[Nsec].flags |= ID_PHOTOM_PASS_0;
+          if (havePS1) GoodPS1 = TRUE;
+          break;
+        case 1:
+          secfilt[Nsec].flags |= ID_PHOTOM_PASS_1;
+          if (havePS1) GoodPS1 = TRUE;
+          break;
+        case 2:
+          secfilt[Nsec].flags |= ID_PHOTOM_PASS_2;
+          if (havePS1) GoodPS1 = TRUE;
+          break;
+        case 3:
+          secfilt[Nsec].flags |= ID_PHOTOM_PASS_3;
+          break;
+        case 4:
+          secfilt[Nsec].flags |= ID_PHOTOM_PASS_4;
+          break;
+      }
+      if (haveSynth) {
+        secfilt[Nsec].flags |= ID_SECF_USE_SYNTH;
+      } 
+      if (havePS1) {
+        secfilt[Nsec].flags |= ID_SECF_HAS_PS1;
+      } 
+      if (haveUbercal) {
+        secfilt[Nsec].flags |= ID_SECF_USE_UBERCAL;
+      } 
+    }
+  }
+
+  if (isSetMrelFinal) {
+    if (pass == 0) {
+      DVOAverageFlags flagBits = ID_OBJ_EXT | ID_OBJ_EXT_ALT | ID_OBJ_GOOD | ID_OBJ_GOOD_ALT;
+      // we attempt to set a few flags here; reset those bits before trying:
+      average[0].flags &= ~flagBits;
+    }
+
+    if (NextPS1 && (NextPS1 > NpsfPS1)) {
+      average[0].flags |= ID_OBJ_EXT;
+    }
+    if (GoodPS1) {
+      average[0].flags |= ID_OBJ_GOOD;
+    }
+    if (Galaxy2MASS) {
+      average[0].flags |= ID_OBJ_EXT_ALT;
+    }
+    if (Good2MASS) {
+      average[0].flags |= ID_OBJ_GOOD_ALT;
+    }
+  }
+  return (TRUE);
+}
+
 // only apply Stack operation on setMrelFinal in first pass 
 // this function has 3 goals, not to be confused: