Changeset 36441 for trunk/psphot

Timestamp:

Jan 14, 2014, 2:28:47 PM (13 years ago)

Author:

bills

Message:

Iniital implementation of the full force stages to the pipeline.
Maginitude and galactic coordinate limits on extended source fits

Location:

trunk

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (2 props)
• psphot/src/Makefile.am (modified) (3 diffs)
• psphot/src/psmakecff.c (modified) (8 diffs)
• psphot/src/psphot.h (modified) (1 diff)
• psphot/src/psphotChooseAnalysisOptions.c (modified) (19 diffs)
• psphot/src/psphotFullForceSummary.c (copied) (copied from branches/bills_branches/bills_201312/psphot/src/psphotFullForceSummary.c )
• psphot/src/psphotFullForceSummaryReadout.c (copied) (copied from branches/bills_branches/bills_201312/psphot/src/psphotFullForceSummaryReadout.c )
• psphot/src/psphotGalaxyShape.c (modified) (3 diffs)
• psphot/src/psphotStackImageLoop.c (modified) (1 prop)

trunk/psphot/src

@@ -26 +26 @@
 psphotFullForce
 psmakecff
+psphotFullForceSummary

@@ -26 +26 @@
 psphotFullForce
 psmakecff
+psphotFullForceSummary

trunk/psphot/src/Makefile.am

@@ -25 +25 @@
 libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 
-bin_PROGRAMS = psphot psphotForced psphotFullForce psphotMinimal psphotMakePSF psphotStack psphotModelTest psmakecff
+bin_PROGRAMS = psphot psphotForced psphotFullForce psphotFullForceSummary psphotMinimal psphotMakePSF psphotStack psphotModelTest psmakecff
 # bin_PROGRAMS = psphotPetrosianStudy psphotTest psphotMomentsStudy 
 
@@ -39 +39 @@
 psphotFullForce_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 psphotFullForce_LDADD = libpsphot.la
+
+psphotFullForceSummary_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
+psphotFullForceSummary_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
+psphotFullForceSummary_LDADD = libpsphot.la
 
 psphotMinimal_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
@@ -99 +103 @@
         psphotMosaicChip.c         \
         psphotCleanup.c
+
+# combine full force results from several inputs
+psphotFullForceSummary_SOURCES = \
+        psphotFullForceSummary.c \
+        psphotFullForceSummaryReadout.c
 
 # forced photometry of specified positions given a specified psf

trunk/psphot/src/psmakecff.c

@@ -2 +2 @@
 # include <config.h>
 # endif
+
+// psmakecff : A program to make read a cmf file and write a cff file
+// The real work is done in psModules.
 
 #include <stdio.h>
@@ -8 +11 @@
 #include "psphot.h"
 
+// For simplicilty, this program's (simple) functions are all contained in this file.
 static pmConfig* psmakecffArguments(int, char**);
 static bool psmakecffParseCamera(pmConfig *);
@@ -25 +29 @@
 //    psphotVersionPrint();
 
-    // load input data (config and images (signal, noise, mask)
     if (!psmakecffParseCamera (config)) {
         psErrorStackPrint(stderr, "Error setting up the camera\n");
@@ -31 +34 @@
     }
 
-    // call psphot for each readout
     if (!psmakecffImageLoop (config)) {
         psErrorStackPrint(stderr, "Error in the psphot image loop\n");
@@ -37 +39 @@
     }
 
-//    psLogMsg ("psphot", PS_LOG_WARN, "complete psphot run: %f sec\n", psTimerMark ("complete"));
-
-//    psErrorCode exit_status = psphotGetExitStatus();
-//    psphotCleanup (config);
     exit (0);
 }
@@ -73 +71 @@
     return config;
 }
+
 static bool psmakecffParseCamera(pmConfig *config) {
     bool status = false;
@@ -127 +126 @@
     psAssert (recipe, "missing recipe?");
 
-//    psImageMaskType maskTest = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
 
     // for psphot, we force data to be read at the chip level
     while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
-        psLogMsg ("psmakecmf", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+        psLogMsg ("psmakecff", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
         if (! chip->process || ! chip->file_exists) { continue; }
 
-#ifdef notdef
-        pmFPAfileActivate (config->files, false, NULL);
-        pmFPAfileActivate (config->files, true, "PSPHOT.LOAD");
-#endif
         if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Chip in psmakecff.");
 
         // there is now only a single chip (multiple readouts?). loop over it and process
         while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
-            psLogMsg ("psmakecmf", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+            psLogMsg ("psmakecff", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
 
             // process each of the readouts
             while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) {
-                psLogMsg ("psmakecmf", 6, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+                psLogMsg ("psmakecff", 6, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
                 if (! readout->data_exists) { continue; }
 
                 pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell);
                 if (hdu && hdu != lastHDU) {
+                    // XXX: probably should do this
                     // psphotVersionHeaderFull(hdu->header);
                     lastHDU = hdu;
@@ -157 +152 @@
 
         }
-        // Defer output and closing of files until we've (possibly) done the NFrames analysis below
-//        pmFPAfileActivate (config->files, false, NULL);
-
         if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed pmFPAfileIOChecks for Chip in psmakecff.");
     }
     if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed pmFPAfileIOChecks for FPA in psphot.");
 
-    // fail if we failed to handle an error
+    // fail if we encountered an unhandled error
     if (psErrorCodeLast() != PS_ERR_NONE) psAbort ("failed to handle an error!");
 

trunk/psphot/src/psphot.h

@@ -361 +361 @@
 bool psphotModelTestReadout(pmConfig *config, const pmFPAview *view, const char *filerule);
 
+bool psphotFullForceSummaryReadout (pmConfig * config, const pmFPAview *view);
+
 int psphotFileruleCount(const pmConfig *config, const char *filerule);
 bool psphotFileruleCountSet(const pmConfig *config, const char *filerule, int num);

trunk/psphot/src/psphotChooseAnalysisOptions.c

@@ -60 +60 @@
 }
 
-bool GetGalacticCoords (psSphere *ptGal, psSphere *ptSky, psSphereRot *toGal, pmChip *chip, float xPos, float yPos);
+static bool GetGalacticCoords (psSphere *ptGal, psSphere *ptSky, psSphereRot *toGal, pmChip *chip, float xPos, float yPos);
 
 // this function use an internal flag to mark sources which have already been measured
@@ -121 +121 @@
     psSphereRot *toGal = NULL;
     float GAL_LIMIT = 0.0;
+    float GAL_LIMIT_BULGE = 0.0;
+    float GAL_LIMIT_SIGMA2 = 0.0;
     bool useGAL_LIMIT = psMetadataLookupBool(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.USE");
     assert (status);
     if (useGAL_LIMIT) {
+        toGal = psSphereRotICRSToGalactic();
         GAL_LIMIT = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT");
-        toGal = psSphereRotICRSToGalactic();
+        assert (status);
+        GAL_LIMIT = DEG_TO_RAD(GAL_LIMIT);
+        GAL_LIMIT_BULGE = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.BULGE");
+        assert (status);
+        GAL_LIMIT_BULGE = DEG_TO_RAD(GAL_LIMIT_BULGE);
+        float GAL_LIMIT_SIGMA = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.BULGE.SIGMA");
+        assert (status);
+        assert(GAL_LIMIT_SIGMA > 0);
+        GAL_LIMIT_SIGMA = DEG_TO_RAD(GAL_LIMIT_SIGMA);
+        GAL_LIMIT_SIGMA2 = GAL_LIMIT_SIGMA * GAL_LIMIT_SIGMA;
     }
 
@@ -135 +147 @@
     }
 
+    float petroFluxLim = NAN;
     float extFitFluxLim = NAN;
 
@@ -140 +153 @@
         float extFitMagLimDefault = NAN;
         float extFitMagLim = NAN;
+        float petroMagLimDefault = NAN;
+        float petroMagLim = NAN;
 
         // match to the given filter
@@ -156 +171 @@
             // find a matching filter or default to 'any'
             if (!strcasecmp (thisFilter, "any")) {
-                psString extFitMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT");
-                psAssert(extFitMagLimStr, "missing MAG.LIMIT");
+                psString petroMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.PETRO");
+                psAssert(petroMagLimStr, "missing MAG.LIMIT.PETRO");
+                petroMagLimDefault = atof (petroMagLimStr);
+
+                psString extFitMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.EXTFIT");
+                psAssert(extFitMagLimStr, "missing MAG.LIMIT.EXTFIT");
                 extFitMagLimDefault = atof (extFitMagLimStr);
             }
 
-            // find a matching filter or default to 'any'
             if (!strcasecmp (thisFilter, filterID)) {
-                psString extFitMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT");
-                psAssert(extFitMagLimStr, "missing MAG.LIMIT");
+                psString petroMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.PETRO");
+                psAssert(petroMagLimStr, "missing MAG.LIMIT.PETRO");
+                petroMagLim = atof (petroMagLimStr);
+
+                psString extFitMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.EXTFIT");
+                psAssert(extFitMagLimStr, "missing MAG.LIMIT.EXTFIT");
                 extFitMagLim = atof (extFitMagLimStr);
                 break;
             }
         }
+        psFree(iter);
+        if (!isfinite (petroMagLim)) petroMagLim = petroMagLimDefault;
         if (!isfinite (extFitMagLim)) extFitMagLim = extFitMagLimDefault;
 
@@ -185 +209 @@
         psAssert (status, "missing FPA.ZP?");
 
+        petroFluxLim = exptime * pow (10.0, 0.4*(zeropt - petroMagLim));
         extFitFluxLim = exptime * pow (10.0, 0.4*(zeropt - extFitMagLim));
     }
@@ -225 +250 @@
             psSphere ptGal, ptSky;
             GetGalacticCoords (&ptGal, &ptSky, toGal, chip, source->peak->xf, source->peak->yf);
-            if (fabs(ptGal.d) < GAL_LIMIT) continue;
+            float b = ptGal.r;
+            float limit = GAL_LIMIT + GAL_LIMIT_BULGE * exp(-0.5*(b*b/GAL_LIMIT_SIGMA2));
+            if (fabs(ptGal.d) < limit) continue;
             // include an exception for low density skycells below the limit?
         }
 
         // for petro and extFit, we will either use the mag limits or the S/N
+        if (doPetrosian) {
+            if (isfinite(petroFluxLim)) {
+                if (source->moments->KronFlux > petroFluxLim) {
+                    source->tmpFlags |= PM_SOURCE_TMPF_PETRO;
+                }
+            } else if (source->moments->KronFlux > SN_LIM * source->moments->KronFluxErr) {
+                source->tmpFlags |= PM_SOURCE_TMPF_PETRO;
+            }
+        }
         if (isfinite(extFitFluxLim)) {
             if (source->moments->KronFlux > extFitFluxLim) {
                 source->tmpFlags |= PM_SOURCE_TMPF_EXT_FIT;
-                if (doPetrosian) source->tmpFlags |= PM_SOURCE_TMPF_PETRO;
             }
         } else {
             if (source->moments->KronFlux > SN_LIM * source->moments->KronFluxErr) {
                 source->tmpFlags |= PM_SOURCE_TMPF_EXT_FIT;
-                if (doPetrosian) source->tmpFlags |= PM_SOURCE_TMPF_PETRO;
             }
         }
@@ -244 +278 @@
 
     psLogMsg ("psphot.options", PS_LOG_WARN, "choose analysis options for %ld sources: %f sec\n", sources->n, psTimerMark ("psphot.options"));
+
+    psFree(toGal);
 
     return true;
@@ -287 +323 @@
     psSphereRot *toGal = NULL;
     float GAL_LIMIT = 0.0;
+    float GAL_LIMIT_BULGE = 0.0;
+    float GAL_LIMIT_SIGMA2 = 0.0;
     bool useGAL_LIMIT = psMetadataLookupBool(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.USE");
     assert (status);
     if (useGAL_LIMIT) {
+        toGal = psSphereRotICRSToGalactic();
         GAL_LIMIT = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT");
-        toGal = psSphereRotICRSToGalactic();
+        assert (status);
+        GAL_LIMIT = DEG_TO_RAD(GAL_LIMIT);
+        GAL_LIMIT_BULGE = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.BULGE");
+        assert (status);
+        GAL_LIMIT_BULGE = DEG_TO_RAD(GAL_LIMIT_BULGE);
+        float GAL_LIMIT_SIGMA = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.BULGE.SIGMA");
+        assert (status);
+        assert(GAL_LIMIT_SIGMA > 0);
+        GAL_LIMIT_SIGMA = DEG_TO_RAD(GAL_LIMIT_SIGMA);
+        GAL_LIMIT_SIGMA2 = GAL_LIMIT_SIGMA * GAL_LIMIT_SIGMA;
     }
 
@@ -365 +413 @@
     // find extFitFluxLim->data.F32[i] for i == image number
     psVector *extFitFluxLim = NULL;
+    psVector *petroFluxLim = NULL;
     if (magLimits) {
         extFitFluxLim = psVectorAlloc (inputFilters->n, PS_TYPE_F32);
+        petroFluxLim = psVectorAlloc (inputFilters->n, PS_TYPE_F32);
         psVectorInit (extFitFluxLim, NAN);
+        psVectorInit (petroFluxLim, NAN);
 
         float extFitMagLimDefault = NAN;
+        float petroMagLimDefault = NAN;
 
         // match mag limits (flux limits) to the filters
@@ -383 +435 @@
             // save the default value to assign to unset filters
             if (!strcasecmp (thisFilter, "any")) {
-                psString magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT");
-                psAssert(magString, "missing MAG.LIMIT");
+                psString magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.EXTFIT");
+                psAssert(magString, "missing MAG.LIMIT.EXTFIT");
                 extFitMagLimDefault = atof (magString);
+
+                magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.PETRO");
+                psAssert(magString, "missing MAG.LIMIT.PETRO");
+                petroMagLimDefault = atof (magString);
                 continue;
             }
@@ -393 +449 @@
                 if (i == chisqNum) continue;
                 if (!strcasecmp (thisFilter, inputFilters->data[i])) {
-                    psString magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT");
-                    psAssert(magString, "missing MAG.LIMIT");
+                    psString magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.PETRO");
+                    psAssert(magString, "missing MAG.LIMIT.PETRO");
                     float magvalue = atof (magString);
+                    petroFluxLim->data.F32[i] = exptime->data.F32[i] * pow (10.0, 0.4*(zeropt->data.F32[i] - magvalue));
+
+                    magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT.EXTFIT");
+                    psAssert(magString, "missing MAG.LIMIT.EXTFIT");
+                    magvalue = atof (magString);
                     extFitFluxLim->data.F32[i] = exptime->data.F32[i] * pow (10.0, 0.4*(zeropt->data.F32[i] - magvalue));
                     break;
@@ -401 +462 @@
             }
         }
+        psFree(iter);
 
         for (int i = 0; i < num; i++) {
             if (i == chisqNum) continue;
-            if (isfinite(extFitFluxLim->data.F32[i])) continue;
-            extFitFluxLim->data.F32[i] = exptime->data.F32[i] * pow (10.0, 0.4*(zeropt->data.F32[i] - extFitMagLimDefault));
+            if (!isfinite(petroFluxLim->data.F32[i])) {
+                petroFluxLim->data.F32[i] = exptime->data.F32[i] * pow (10.0, 0.4*(zeropt->data.F32[i] - petroMagLimDefault));
+            }
+            if (!isfinite(extFitFluxLim->data.F32[i])) {
+                extFitFluxLim->data.F32[i] = exptime->data.F32[i] * pow (10.0, 0.4*(zeropt->data.F32[i] - extFitMagLimDefault));
+            }
         }
     }
@@ -420 +486 @@
         bool doObjectRadial = false;
         bool doObjectExtFit = false;
-        for (int j = 0; !doObjectExtFit && !doObjectRadial && (j < object->sources->n); j++) {
+        bool doObjectPetrosian = false;
+        for (int j = 0; !doObjectExtFit && !doObjectRadial && !doObjectPetrosian && (j < object->sources->n); j++) {
 
             pmSource *source = object->sources->data[j];
@@ -451 +518 @@
             if (extFitAll) {
                 doObjectExtFit = true;
+                doObjectPetrosian = doPetrosian;
                 continue;
             }
@@ -465 +533 @@
                 psSphere ptGal, ptSky;
                 GetGalacticCoords (&ptGal, &ptSky, toGal, chips->data[imageID], source->peak->xf, source->peak->yf);
-                if (fabs(ptGal.d) < GAL_LIMIT) continue;
-                // include an exception for low density skycells below the limit?
-            }
-
-            float fluxLim = extFitFluxLim ? extFitFluxLim->data.F32[imageID] : NAN;
+                float b = ptGal.r;
+                float limit = GAL_LIMIT + GAL_LIMIT_BULGE * exp(-0.5*(b*b/GAL_LIMIT_SIGMA2));
+                if (fabs(ptGal.d) < limit) continue;
+            }
+
+            float fluxLim = NAN;
+            // for petro and extFit, we will either use the mag limits or the S/N
+            if (doPetrosian) {
+                fluxLim = petroFluxLim ? petroFluxLim->data.F32[imageID] : NAN;
+                if (isfinite(fluxLim)) {
+                    if (source->moments->KronFlux > extFitFluxLim->data.F32[imageID]) {
+                        doObjectPetrosian = true;
+                    }
+                } else {
+                    if (source->moments->KronFlux > SN_LIM * source->moments->KronFluxErr) {
+                        doObjectPetrosian = true;
+                    }
+                }
+            }
+            fluxLim = extFitFluxLim ? extFitFluxLim->data.F32[imageID] : NAN;
             // for petro and extFit, we will either use the mag limits or the S/N
             if (isfinite(fluxLim)) {
@@ -503 +586 @@
             if (source->modelPSF == NULL) continue;
                 
-            // Do the fits if the recipe requests we do extended source fits to everything
             if (doObjectExtFit) {
                 source->tmpFlags |= PM_SOURCE_TMPF_EXT_FIT;
-                if (doPetrosian) source->tmpFlags |= PM_SOURCE_TMPF_PETRO;
-            }
+            }
+            if (doObjectPetrosian) {
+                source->tmpFlags |= PM_SOURCE_TMPF_PETRO;
+            }
 
             // Do the fits if the recipe requests we do extended source fits to everything
@@ -522 +606 @@
     psLogMsg ("psphot.options", PS_LOG_WARN, "choose analysis options for %ld objects: %f sec\n", objects->n, psTimerMark ("psphot.options"));
 
+    psFree(toGal);
+
     return true;
 }
 
-bool GetGalacticCoords (psSphere *ptGal, psSphere *ptSky, psSphereRot *toGal, pmChip *chip, float xPos, float yPos) {
+static bool GetGalacticCoords (psSphere *ptGal, psSphere *ptSky, psSphereRot *toGal, pmChip *chip, float xPos, float yPos) {
 
     pmFPA *fpa = chip->parent;
@@ -543 +629 @@
     psSphereRotApply (ptGal, toGal, ptSky);
 
+    // psSphereRotApply insures that 0 < r < 2PI. We want -PI < b <= PI
+    if (ptGal->r > M_PI) {
+        ptGal->r -= 2.0 * M_PI;
+    }
+
     return true;
 

trunk/psphot/src/psphotGalaxyShape.c

@@ -283 +283 @@
     }
 
+#ifdef SAVE_BEST_MODEL
+    // Save model with smallest chisq
     if (isfinite(chisqBest)) {
+#else 
+    // Save model with nominal parameters
+    {
+        fRmajorBest = 1;
+        fRminorBest = 1;
+#endif
         // now save the best fitting model as the source's extended model
         psEllipseAxes testAxes = guessAxes;
@@ -293 +301 @@
         psphotGalaxyShapeSource (pcm, source, maskVal, psfSize, false);
 
+        // Replace modelEXT with this model
+        // XXX: only do this if the model is good
         psFree (source->modelEXT);
 
@@ -299 +309 @@
         source->mode |= PM_SOURCE_MODE_EXTMODEL;
         source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
-
-        // adjust the window so the subtraction covers the faint wings
-        // psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal);
 
         // cache the model flux

@@ -26 +26 @@
 psphotFullForce
 psmakecff
+psphotFullForceSummary