Changeset 25204

Timestamp:

Aug 26, 2009, 7:49:18 AM (17 years ago)

Author:

eugene

Message:

adding errors to the stats

File:

• branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c (modified) (3 diffs)

branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c

@@ -18 +18 @@
 
     psVector *fluxSum     = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
+    psVector *fluxSumErr2 = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
     psVector *refRadius   = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
     psVector *petRatio    = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
@@ -27 +28 @@
     bool manyPetro = false;
     bool above = true;
+    float Asum = 0.0;
     float Fsum = 0.0;
-    float Asum = 0.0;
+    float dFsum2 = 0.0;
 
     int nOut = 0;
@@ -36 +38 @@
 
         float Area = area->data.F32[i];
+        Asum += Area;
         Fsum += binSB->data.F32[i] * Area;
-        Asum += Area;
+        dFsum2 += PS_SQR(binSBstdev->data.F32[i] * Area);
 
+        float areaInner = 0.5 * Area;
         float fluxInner = 0.5 * Area * binSB->data.F32[i];
-        float areaInner = 0.5 * Area;
+        float fluxInnerErr2 = PS_SQR(binSBstdev->data.F32[i] * 0.5 * Area);
         if (nOut > 0) {
+            areaInner += areaSum->data.F32[nOut-1];
             fluxInner += fluxSum->data.F32[nOut-1];
-            areaInner += areaSum->data.F32[nOut-1];
+            fluxInnerErr2 = PS_SQR(fluxSumErr2->data.F32[nOut-1] * areaSum->data.F32[nOut-1]);
         }
 
+        // ratio = binSB / meanSB
+        // meanSB = flux / area
+        // flux = sum(binSB(i) * area(i)
+        // fluxErr^2 = sum(binSBerr(i)^2 area(i)^2)
+        // meanSBerr^2 = fluxErr^2 / area^2
+        // (ratioErr/ratio)^2 = (binSBerr/binSB)^2 + (meanSBerr/meanSB)^2
+
         psVectorAppend(meanSB, (fluxInner / areaInner));
-        psVectorAppend(petRatio, binSB->data.F32[i] / meanSB->data.F32[nOut]);
-        psVectorAppend(petRatioErr, binSBstdev->data.F32[i] / meanSB->data.F32[nOut]);
+
+        float ratio = binSB->data.F32[i] / meanSB->data.F32[nOut];
+        psVectorAppend(petRatio, ratio);
+
+        float meanSBerr = sqrt(fluxInnerErr2) / areaInner;
+        float ratioErr = ratio * sqrt(PS_SQR(binSBstdev->data.F32[i]/binSB->data.F32[i]) + PS_SQR(meanSBerr/meanSB->data.F32[nOut]));
+
+        psVectorAppend(petRatioErr, ratioErr);
+
+        psVectorAppend(areaSum, Asum);
         psVectorAppend(fluxSum, Fsum);
-        psVectorAppend(areaSum, Asum);
+        psVectorAppend(fluxSumErr2, dFsum2);
         psVectorAppend(refRadius, binRad->data.F32[i]);