Changeset 9834 for trunk/psphot/src/psphotEnsemblePSF.c

Timestamp:

Nov 3, 2006, 5:05:49 AM (20 years ago)

Author:

eugene

Message:

fixed formula for normalization error

File:

• trunk/psphot/src/psphotEnsemblePSF.c (modified) (4 diffs)

trunk/psphot/src/psphotEnsemblePSF.c

@@ -141 +141 @@
     psSparse *sparse = psSparseAlloc (models->n, 100);
     psVector *weight = psVectorAlloc (models->n, PS_TYPE_F32);
+    psVector *errors = psVectorAlloc (models->n, PS_TYPE_F32);
 
     for (int i = 0; i < models->n; i++) {
@@ -150 +151 @@
         r = pmSourceCrossProduct (Mi, Mi, CONSTANT_PHOTOMETRIC_WEIGHTS);
         weight->data.F32[i] = r;
+
+        // the formal error depends on the weighting scheme
+        if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
+            float var = pmSourceCrossProduct (Mi, Mi, false);
+            errors->data.F32[i] = 1.0 / sqrt(var);
+        } else {
+            errors->data.F32[i] = 1.0 / sqrt(r);
+        }
 
         psSparseMatrixElement (sparse, i, i, 1.0);
@@ -200 +209 @@
         }
         Fi->modelPSF->params->data.F32[PM_PAR_I0] = norm->data.F32[i];
-        Fi->modelPSF->dparams->data.F32[PM_PAR_I0] = sqrt(sqrt(2) * norm->data.F32[i] / (sparse->Bfj->data.F32[i] * weight->data.F32[i]));
+        // XXX I was using this before, but it is not strictly accurate
         // XXX EAM : this factor of sqrt(2) makes the errors consistent, but I don't understand it
+        // Fi->modelPSF->dparams->data.F32[PM_PAR_I0] = sqrt(sqrt(2) * norm->data.F32[i] / (sparse->Bfj->data.F32[i] * weight->data.F32[i]));
+        Fi->modelPSF->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i];
 
         // subtract object
@@ -228 +239 @@
     psFree (norm);
     psFree (weight);
+    psFree (errors);
 
     psLogMsg ("psphot.emsemble", 3, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot"));