Changeset 40231

Timestamp:

Nov 30, 2017, 8:45:51 PM (9 years ago)

Author:

eugene

Message:

need to set Mkron for average chips or we cannot measure the stack zp.aper

Location:

branches/eam_branches/ohana.20170822/src/relphot

Files:

• include/relphot.h (modified) (1 diff)
• src/StarOps.c (modified) (1 diff)
• src/relphot_images.c (modified) (1 diff)
• src/setMrelCatalog.c (modified) (6 diffs)

branches/eam_branches/ohana.20170822/src/relphot/include/relphot.h

@@ -406 +406 @@
 void          findImages          PROTO((Catalog *catalog, int Ncatalog, int doImageList));
 int           findMosaics         PROTO((Catalog *catalog, int Ncatalog, int doMosaicList));
+
+int dumpMags (FILE *fout, Catalog *catalog, int Ncatalog);
 
 void makeMosaics (Image *image, off_t Nimage, int mergeMcal);

branches/eam_branches/ohana.20170822/src/relphot/src/StarOps.c

@@ -542 +542 @@
   }
   fclose (ftest);
+  return TRUE;
+}
+
+int dumpMags (FILE *fout, Catalog *catalog, int Ncatalog) {
+
+  int i, n;
+  off_t j;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  for (i = 0; i < Ncatalog; i++) {
+    for (j = 0; (j < catalog[i].Naverage) && (j < 2); j++) {
+      fprintf (fout, "%08x %08x %10.6f %10.6f : \n", catalog[i].averageT[j].catID, catalog[i].averageT[j].objID, catalog[i].averageT[j].R, catalog[i].averageT[j].D); 
+      for (n = 0; n < 5; n++) {
+        fprintf (fout, "secf  %5d | %6.3f %6.3f | %6.3f %6.3f | %6.3f %6.3f\n", n, catalog[i].secfilt[j*Nsecfilt + n].MpsfChp, catalog[i].secfilt[j*Nsecfilt + n].MkronChp, catalog[i].secfilt[j*Nsecfilt + n].MpsfStk, catalog[i].secfilt[j*Nsecfilt + n].MkronStk, catalog[i].secfilt[j*Nsecfilt + n].MpsfWrp, catalog[i].secfilt[j*Nsecfilt + n].MkronWrp);
+      }
+      int m = catalog[i].averageT[j].measureOffset;
+      for (n = 0; n < catalog[i].averageT[j].Nmeasure; n++) {
+        fprintf (fout, "meas %5d %5d | %6.3f %6.3f | %6.3f\n", m+n, catalog[i].measureT[m+n].photcode, catalog[i].measureT[m+n].M, catalog[i].measureT[m+n].Mkron, catalog[i].measureT[m+n].Mcal); 
+      }
+    }
+  }
   return TRUE;
 }

branches/eam_branches/ohana.20170822/src/relphot/src/relphot_images.c

@@ -108 +108 @@
         plot_scatter (catalog, Ncatalog, flatcorr); 
       }
+      if (VERBOSE2) dumpMags (stderr, catalog, Ncatalog);
+
       setMrel  (catalog, Ncatalog, flatcorr); // threaded
+      if (VERBOSE2) dumpMags (stderr, catalog, Ncatalog);
+
       if (PLOTSTUFF) {
         plot_scatter (catalog, Ncatalog, flatcorr); 
       }
       setMcal  (catalog, FALSE, flatcorr);
+      if (VERBOSE2) dumpMags (stderr, catalog, Ncatalog);
+
       setMmos  (catalog, FALSE, flatcorr);
       setMgrid (catalog, flatcorr);

branches/eam_branches/ohana.20170822/src/relphot/src/setMrelCatalog.c

@@ -267 +267 @@
     }
 
-    float Mkron, dMkron;
+    // NOTE: we need to calculate the averge chip Kron on each pass to be able to calibrate the stacks
+    float Mkron = PhotCatTiny (&measureT[k], MAG_CLASS_KRON);
+    float dMkron;
     if (isSetMrelFinal) {
-      Mkron = PhotCat (&measure[k], MAG_CLASS_KRON);
       dMkron = MAX (hypot(measure[k].dMkron, code->photomErrSys), MIN_ERROR);
-      if (CHECK_VALID_MAG(Mkron, dMkron)) {
-        int Nkron = results->kronData[Nsec].Nlist;
-        results->kronData[Nsec].flxlist[Nkron] = Mkron - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel
-        results->kronData[Nsec].errlist[Nkron] = dMkron;
-        results->kronData[Nsec].wgtlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0;
-        results->kronData[Nsec].ranking[Nkron] = measureRank[k];
-        results->kronData[Nsec].measSeq[Nkron] = k;
-        results->kronData[Nsec].msklist[Nkron] = 0;
-        results->kronData[Nsec].Nlist ++;
-      }
+    } else {
+      dMkron = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR);
+    }
+    if (CHECK_VALID_MAG(Mkron, dMkron)) {
+      int Nkron = results->kronData[Nsec].Nlist;
+      results->kronData[Nsec].flxlist[Nkron] = Mkron - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel
+      results->kronData[Nsec].errlist[Nkron] = dMkron;
+      results->kronData[Nsec].wgtlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0;
+      results->kronData[Nsec].ranking[Nkron] = measureRank ? measureRank[k] : 0;
+      results->kronData[Nsec].measSeq[Nkron] = k;
+      results->kronData[Nsec].msklist[Nkron] = 0;
+      results->kronData[Nsec].Nlist ++;
     }
 
@@ -385 +388 @@
     // if too few valid measurements meet the minimum criteria, go to the next entry
     StatType *psfstats = &results->psfstats;
-    int Nranking = MAG_STATS_BY_RANKING (&results->psfData[Nsec], psfstats);
-    if (Nranking < Nminmeas) { 
+    int NrankingPSF = MAG_STATS_BY_RANKING (&results->psfData[Nsec], psfstats);
+    if (NrankingPSF < Nminmeas) { 
       secfilt[Nsec].flags |= ID_OBJ_FEW;
     } else {
@@ -393 +396 @@
       secfilt[Nsec].Mchisq   = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
     }
-    int minRankPSF = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10;
 
     // when running -averages, we have no information about the images, so we cannot set this
@@ -400 +402 @@
     }
 
+    StatType *kronstats = &results->kronstats;
+    int NrankingKRON = MAG_STATS_BY_RANKING (&results->kronData[Nsec], kronstats);
+    if (NrankingKRON) {
+      secfilt[Nsec].MkronChp  = kronstats->mean; 
+      secfilt[Nsec].dMkronChp = kronstats->error; 
+      secfilt[Nsec].sMkronChp = kronstats->sigma; 
+      secfilt[Nsec].NusedKron = NrankingKRON;
+    }
+
     if (isSetMrelFinal) {
       if ((average[0].objID == TEST_OBJ_ID) && (average[0].catID == TEST_CAT_ID)) {
@@ -405 +416 @@
       }
 
-      // mark the measurements matching this ranking 
-      markMeasureByRanking (&results->psfData[Nsec], measure, minRankPSF, ID_MEAS_PHOTOM_PSF, ID_MEAS_MASKED_PSF);
-
-      if (Nranking) {
+      if (NrankingPSF) {
         secfilt[Nsec].sMpsfChp = psfstats->sigma;
         secfilt[Nsec].Nused    = psfstats->Nmeas;
@@ -420 +428 @@
       MAX_NOT_NAN (psfQfPerfMax, secfilt[Nsec].psfQfPerfMax);
 
+      // mark the measurements matching this ranking 
+      int minRankPSF = (NrankingPSF > 0) ? results->psfData[Nsec].ranking[0] : 10;
+      markMeasureByRanking (&results->psfData[Nsec], measure, minRankPSF, ID_MEAS_PHOTOM_PSF, ID_MEAS_MASKED_PSF);
+
+      int minRankKron = (NrankingKRON > 0) ? results->kronData[Nsec].ranking[0] : 10;
+      markMeasureByRanking (&results->kronData[Nsec], measure, minRankKron, ID_MEAS_PHOTOM_KRON, ID_MEAS_MASKED_KRON);
+
       StatType *apstats = &results->apstats;
-      Nranking = MAG_STATS_BY_RANKING (&results->aperData[Nsec], apstats);
-      if (Nranking) {
+      int NrankingAPER = MAG_STATS_BY_RANKING (&results->aperData[Nsec], apstats);
+      if (NrankingAPER) {
         secfilt[Nsec].MapChp  = apstats->mean; 
         secfilt[Nsec].dMapChp = apstats->error; 
         secfilt[Nsec].sMapChp = apstats->sigma; 
-        secfilt[Nsec].NusedAp = Nranking;
-      }
-      int minRankAper = (Nranking > 0) ? results->aperData[Nsec].ranking[0] : 10;
+        secfilt[Nsec].NusedAp = NrankingAPER;
+      }
+      int minRankAper = (NrankingAPER > 0) ? results->aperData[Nsec].ranking[0] : 10;
       markMeasureByRanking (&results->aperData[Nsec], measure, minRankAper, ID_MEAS_PHOTOM_APER, ID_MEAS_MASKED_APER);
-
-      StatType *kronstats = &results->kronstats;
-      Nranking = MAG_STATS_BY_RANKING (&results->kronData[Nsec], kronstats);
-      if (Nranking) {
-        secfilt[Nsec].MkronChp  = kronstats->mean; 
-        secfilt[Nsec].dMkronChp = kronstats->error; 
-        secfilt[Nsec].sMkronChp = kronstats->sigma; 
-        secfilt[Nsec].NusedKron = Nranking;
-      }
-      int minRankKron = (Nranking > 0) ? results->kronData[Nsec].ranking[0] : 10;
-      markMeasureByRanking (&results->kronData[Nsec], measure, minRankKron, ID_MEAS_PHOTOM_KRON, ID_MEAS_MASKED_KRON);
 
       // does this object appear extended in > 50% of measurements?