Changeset 36747 for trunk/psphot/src/psphotGalaxyShape.c

Timestamp:

May 13, 2014, 2:40:41 PM (12 years ago)

Author:

bills

Message:

change psphotFullForce galaxy shapes analysis to use different radius
trial values based on the signal to nosie of the source as measured
by psphotStack

File:

• trunk/psphot/src/psphotGalaxyShape.c (modified) (7 diffs)

trunk/psphot/src/psphotGalaxyShape.c

@@ -80 +80 @@
     assert (markVal);
 
-    // maskVal is used to test for rejected pixels, and must include markVal
     maskVal |= markVal;
 
-    // bit-mask to mark pixels not used in analysis
     psphotGalaxyShapeOptions *opt = psphotGalaxyShapeOptionsAlloc();
+    opt->Q = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_Q"); 
+    psAssert (status, "missing GALAXY_SHAPES_Q");
+    opt->NSigma = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_NSIGMA"); 
+    psAssert (status, "missing GALAXY_SHAPES_NSIGMA");
+    opt->clampSN = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_CLAMP_SN"); 
+    psAssert (status, "missing GALAXY_SHAPES_NSIGMA");
+
+#ifdef notdef
     opt->fRmajorMin = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MAJOR_MIN"); if (!status) opt->fRmajorMin = 0.5;
     opt->fRmajorMax = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MAJOR_MAX"); if (!status) opt->fRmajorMax = 2.0;
@@ -91 +97 @@
     opt->fRminorMax = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MINOR_MAX"); if (!status) opt->fRminorMax = 2.0;
     opt->fRminorDel = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MINOR_DEL"); if (!status) opt->fRminorDel = 0.1;
+#endif
 
     // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
@@ -238 +245 @@
 bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psphotGalaxyShapeOptions *opt, psImageMaskType maskVal, int psfSize) {
 
+
     for (int iModel = 0 ; iModel < source->modelFits->n; iModel++) {
         pmModel *model = source->modelFits->data[iModel];
@@ -252 +260 @@
 
         // double check that the guess is carried along...
-
-        // I have some source guess (e0, e1, e2)
-        psEllipseAxes guessAxes = pmPSF_ModelToAxes (PAR, modelType);
-
         if (!source->galaxyFits) {
             source->galaxyFits = psArrayAllocEmpty(1);
-            // source->galaxyFits = pmSourceGalaxyFitsAlloc();
         }
         pmSourceGalaxyFits *galaxyFits = pmSourceGalaxyFitsAlloc();
@@ -264 +267 @@
         psFree(galaxyFits);
 
+        // set the options for this source
+        psphotGalaxyShapeOptionsSet(source, galaxyFits, opt);
+
+        // I have some source guess (e0, e1, e2)
+        psEllipseAxes guessAxes = pmPSF_ModelToAxes (PAR, modelType);
+
         float fRmajorBest = NAN;
         float fRminorBest = NAN;
         float chisqBest = NAN;
-        for (float fRmajor = opt->fRmajorMin; fRmajor < opt->fRmajorMax + 0.5*opt->fRmajorDel; fRmajor += opt->fRmajorDel) {
-            for (float fRminor = opt->fRminorMin; fRminor < opt->fRminorMax + 0.5*opt->fRminorDel; fRminor += opt->fRminorDel) {
+        for (float fRmajor = galaxyFits->fRmajorMin; fRmajor < galaxyFits->fRmajorMax + 0.5*galaxyFits->fRmajorDel; fRmajor += galaxyFits->fRmajorDel) {
+            for (float fRminor = galaxyFits->fRminorMin; fRminor < galaxyFits->fRminorMax + 0.5*galaxyFits->fRminorDel; fRminor += galaxyFits->fRminorDel) {
       
                 psEllipseAxes testAxes = guessAxes;
@@ -405 +414 @@
     psMemSetDeallocator(opt, (psFreeFunc) psphotGalaxyShapeOptionsFree);
     
+#ifdef notdef
     opt->fRmajorMin = 0.5;
     opt->fRmajorMax = 2.0;
@@ -411 +421 @@
     opt->fRminorMax = 2.0;
     opt->fRminorDel = 0.1;
+    opt->snMinAllTrials = 20;      // for sources with cff SN < 20 ...
+    opt->reduceTrialsFactor = 2; // cut the number of trials in half (double the deltas)
+#endif
 
     return opt;
 }
 
+void psphotGalaxyShapeOptionsSet(pmSource *source, pmSourceGalaxyFits *galaxyFits, psphotGalaxyShapeOptions *options)
+{
+    // XXX: put these in recipe
+    // doesn't make sense to use too small of a spacing
+    float clampedSN = source->extSN < options->clampSN ? source->extSN : options->clampSN;
+
+    float f_del = options->Q / clampedSN;
+    float f_min = 1 - options->NSigma * f_del;
+    float f_max = 1 + options->NSigma * f_del;
+
+    galaxyFits->fRmajorMin = f_min;
+    galaxyFits->fRmajorMax = f_max;
+    galaxyFits->fRmajorDel = f_del;
+    galaxyFits->fRminorMin = f_min;
+    galaxyFits->fRminorMax = f_max;
+    galaxyFits->fRminorDel = f_del;
+}
+