Changeset 9860 for trunk/Mars/mimage/MHillasExt.cc

Timestamp:

08/13/10 16:45:16 (14 years ago)

Author:

tbretz

Message:

Use the definition of the biased weighted sample variance to calculate the weighted time spread. This includes the weights (like in the third moment) linearly, not squared as before.

File:

• trunk/Mars/mimage/MHillasExt.cc (modified) (4 diffs)

trunk/Mars/mimage/MHillasExt.cc

@@ -154 +154 @@
     Double_t sumy      = 0;
     Double_t sumt      = 0;
-    Double_t sumt2     = 0;
+    Double_t sumw      = 0;
+
+    Double_t sumxt     = 0;
+    Double_t sumyt     = 0;
+
     Double_t sumx2     = 0;
     Double_t sumy2     = 0;
     Double_t sumxy     = 0;
-    Double_t sumxt     = 0;
-    Double_t sumyt     = 0;
-    Double_t sumw      = 0;
-    Double_t sumw2     = 0;
+    Double_t sumt2     = 0;
+
     Double_t sumxw     = 0;
     Double_t sumyw     = 0;
     Double_t sumtw     = 0;
-    Double_t sumtx2    = 0;
-    Double_t sumty2    = 0;
-    Double_t sumtw2    = 0;
-    Double_t sumxw2    = 0;
-    Double_t sumyw2    = 0;
-    Double_t sumtxw2   = 0;
-    Double_t sumtyw2   = 0;
-    Double_t sumt2w2   = 0;
-    Double_t sumx2w2   = 0;
-    Double_t sumy2w2   = 0;
-    Double_t sumxyw2   = 0;
+
+    Double_t sumx2w    = 0;
+    Double_t sumy2w    = 0;
+    Double_t sumt2w    = 0;
+    Double_t sumxyw    = 0;
+    Double_t sumxtw    = 0;
+    Double_t sumytw    = 0;
 
     Double_t sumdx3w   = 0;
@@ -204 +202 @@
         sumx    += x;
         sumy    += y;
-
         sumt    += t;
-        sumt2   += t*t;
+        sumw    += nphot;
+
+        sumxt   += x*t;
+        sumyt   += y*t;
 
         sumx2   += x*x;
         sumy2   += y*y;
         sumxy   += x*y;
-        sumxt   += x*t;
-        sumyt   += y*t;
-
-        sumw    += nphot;
-        sumw2   += nphot*nphot;
+        sumt2   += t*t;
 
         sumxw   += x*nphot;
@@ -221 +217 @@
         sumtw   += t*nphot;
 
-        sumtx2  += t*x*x;
-        sumty2  += t*y*y;
-        sumtw2  += t*nphot*nphot;
-
-        sumxw2  += x*nphot*nphot;
-        sumyw2  += y*nphot*nphot;
-
-        sumtxw2 += t*x*nphot*nphot;
-        sumtyw2 += t*y*nphot*nphot;
-        sumt2w2 += t*t*nphot*nphot;
-        sumx2w2 += x*x*nphot*nphot;
-        sumy2w2 += y*y*nphot*nphot;
-        sumxyw2 += x*y*nphot*nphot;
+        sumx2w  += x*x*nphot;
+        sumy2w  += y*y*nphot;
+        sumt2w  += t*t*nphot;
+
+        sumxyw  += x*y*nphot;
+
+        sumxtw  += x*t*nphot;
+        sumytw  += y*t*nphot;
 
         // --- 3rd moment ---
@@ -286 +277 @@
     const Double_t m = fSlopeLong;
 
-    const Double_t z        = c*sumx    + s*sumy;
-    const Double_t zw       = c*sumxw   + s*sumyw;
-    const Double_t zw2      = c*sumxw2  + s*sumyw2;
-    const Double_t zt       = c*sumxt   + s*sumyt;
-    const Double_t ztw2     = c*sumtxw2 + s*sumtyw2;
-
-    const Double_t sumdt    = sumt   - m*z;
-    const Double_t sumdtw   = sumtw  - m*zw;
-    const Double_t sumdtw2  = sumtw2 - m*zw2;
-
-    const Double_t x2y2     = c*c*sumx2   + s*s*sumy2   + 2*c*s*sumxy;
-    const Double_t x2y2w2   = c*c*sumx2w2 + s*s*sumy2w2 + 2*c*s*sumxyw2;
-
-    const Double_t sumdt2   = sumt2   + m*m*x2y2   - 2*m*zt;
-    const Double_t sumdt2w2 = sumt2w2 + m*m*x2y2w2 - 2*m*ztw2;
-
-    const Double_t meant    = sumt /cnt;
-    const Double_t meandt   = sumdt/cnt;  // meant - m*z/cnt
-
+    const Double_t z        = dx; //c*sumx   + s*sumy;
+    const Double_t zw       =       c*sumxw  + s*sumyw;
+
+    const Double_t zt       = dxt; //c*sumxt  + s*sumyt;
+    const Double_t ztw      =        c*sumxtw + s*sumytw;
+
+    const Double_t x2y2     = dx2; //c*c*sumx2  + s*s*sumy2  + 2*c*s*sumxy;
+    const Double_t x2y2w    =        c*c*sumx2w + s*s*sumy2w + 2*c*s*sumxyw;
+
+    const Double_t sumdt    = sumt  - m*z;
+    const Double_t sumdtw   = sumtw - m*zw;
+
+    const Double_t sumdt2   = sumt2  + m*m*x2y2  - 2*m*zt;
+    const Double_t sumdt2w  = sumt2w + m*m*x2y2w - 2*m*ztw;
+
+    const Double_t meant    = sumt  /cnt;
+    const Double_t meandt   = sumdt /cnt;
     const Double_t meantw   = sumtw /sumw;
     const Double_t meandtw  = sumdtw/sumw;
 
-    fTimeSpread          = TMath::Sqrt(sumt2 /cnt - meant *meant);
-    fSlopeSpread         = TMath::Sqrt(sumdt2/cnt - meandt*meandt);
-
-    fTimeSpreadWeighted  = TMath::Sqrt(sumt2w2  + sumw2*meantw *meantw  - 2*meantw *sumtw2) /sumw;
-    fSlopeSpreadWeighted = TMath::Sqrt(sumdt2w2 + sumw2*meandtw*meandtw - 2*meandtw*sumdtw2)/sumw;
+    fTimeSpread          = TMath::Sqrt(sumt2 /cnt  - meant *meant);
+    fTimeSpreadWeighted  = TMath::Sqrt(sumt2w/sumw - meantw*meantw);
+
+    fSlopeSpread         = TMath::Sqrt(sumdt2 /cnt  - meandt *meandt);
+    fSlopeSpreadWeighted = TMath::Sqrt(sumdt2w/sumw - meandtw*meandtw);
 
     //