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

Timestamp:

08/13/10 11:34:45 (14 years ago)

Author:

tbretz

Message:

Fixed the calculation of the weighted time spread and a bug in the
calculation of the third moments!

File:

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

trunk/Mars/mimage/MHillasExt.cc

@@ -151 +151 @@
     UInt_t cnt = 0;
 
-    Double_t sumx    = 0;
-    Double_t sumy    = 0;
-    Double_t sumt    = 0;
-    Double_t sumw    = 0;
-    Double_t sumxy   = 0;
-    Double_t sumxt   = 0;
-    Double_t sumyt   = 0;
-    Double_t sumxw   = 0;
-    Double_t sumyw   = 0;
-    Double_t sumtw   = 0;
-    Double_t sumx2   = 0;
-    Double_t sumy2   = 0;
-    Double_t sumw2   = 0;
-    Double_t sumt2   = 0;
-    Double_t sumx3   = 0;
-    Double_t sumy3   = 0;
-    Double_t sumxw2  = 0;
-    Double_t sumyw2  = 0;
-    Double_t sumtx2  = 0;
-    Double_t sumty2  = 0;
-    Double_t sumtw2  = 0;
-    Double_t sumx2y  = 0;
-    Double_t sumxy2  = 0;
-    Double_t sumtxw2 = 0;
-    Double_t sumtyw2 = 0;
-    Double_t sumt2w2 = 0;
+    Double_t sumx      = 0;
+    Double_t sumy      = 0;
+    Double_t sumt      = 0;
+    Double_t sumt2     = 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 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 sumdx3w   = 0;
+    Double_t sumdy3w   = 0;
+
+    Double_t sumdx2dyw = 0;
+    Double_t sumdxdy2w = 0;
 
     const UInt_t npix = evt.GetNumPixels();
@@ -203 +208 @@
         sumt2   += t*t;
 
+        sumx2   += x*x;
+        sumy2   += y*y;
+        sumxy   += x*y;
+        sumxt   += x*t;
+        sumyt   += y*t;
+
         sumw    += nphot;
         sumw2   += nphot*nphot;
@@ -210 +221 @@
         sumtw   += t*nphot;
 
-        sumtx2  += x*nphot*nphot;
-        sumty2  += y*nphot*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;
-
-        sumx2 += x*x;
-        sumy2 += y*y;
-        sumxy += x*y;
-        sumxt += x*t;
-        sumyt += y*t;
+        sumx2w2 += x*x*nphot*nphot;
+        sumy2w2 += y*y*nphot*nphot;
+        sumxyw2 += x*y*nphot*nphot;
 
         // --- 3rd moment ---
@@ -228 +239 @@
         const Double_t dy = y - hil.GetMeanY();      // [mm]
 
-        sumx3  += nphot * dx*dx*dx;
-        sumy3  += nphot * dy*dy*dy;
-        sumx2y += nphot * dx*dx*dy;
-        sumx2y += nphot * dx*dy*dy;
+        sumdx3w   += dx*dx*dx*nphot;
+        sumdy3w   += dy*dy*dy*nphot;
+        sumdx2dyw += dx*dx*dy*nphot;
+        sumdxdy2w += dx*dy*dy*nphot;
 
         cnt++;
@@ -253 +264 @@
 
     //
-    // Time spread
-    //
-    fTimeSpread         = TMath::Sqrt(sumt2*cnt - sumt*sumt)/cnt;
-    fTimeSpreadWeighted = TMath::Sqrt(sumt2w2 + sumw2*sumtw*sumtw/sumw/sumw - 2*sumtw2*sumtw/sumw)/sumw;
-
-    //
     // Time slope
     //
@@ -277 +282 @@
 
     //
-    // Calculate time spread around longitudinal time slope
-    //
-    const Double_t sumdt     = sumt    -   fSlopeLong*(c*sumx    + s*sumy);
-    const Double_t sumdtw    = sumtw   -   fSlopeLong*(c*sumxw   + s*sumyw);
-    const Double_t sumdtw2   = sumtw2  -   fSlopeLong*(c*sumxw2  + s*sumyw2);
-    const Double_t sumdt2    = sumt2   - 2*fSlopeLong*(c*sumxt   + s*sumyt)   + fSlopeLong*fSlopeLong * (c*c*sumx2 + s*s*sumy2 + 2*c*s*sumxy);
-    const Double_t sumdt2w2  = sumt2w2 - 2*fSlopeLong*(c*sumtxw2 + s*sumtyw2) + fSlopeLong*fSlopeLong * (c*c*sumx2 + s*s*sumy2 + 2*c*s*sumxy) * sumw2;
-
-    fSlopeSpread         = TMath::Sqrt(sumdt2*cnt - sumdt*sumdt)/cnt;
-    fSlopeSpreadWeighted = TMath::Sqrt(sumdt2w2 + sumw2*sumdtw*sumdtw/sumw/sumw - 2*sumdtw2*sumdtw/sumw)/sumw;
+    // Time spread
+    //
+    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 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;
 
     //
     // Third moments along axes get normalized
     //
-    const Double_t m3l = c*c*c*sumx3 + 3*(s*c*c*sumx2y + c*s*s*sumxy2) + s*s*s*sumy3;
-    const Double_t m3t = c*c*c*sumy3 + 3*(s*s*c*sumx2y - s*c*c*sumxy2) - s*s*s*sumx3;
-
-    fM3Long  = MMath::Sqrt3(m3l/hil.GetSize());      // [mm]
-    fM3Trans = MMath::Sqrt3(m3t/hil.GetSize());      // [mm]
+    const Double_t m3l = c*c*c*sumdx3w + s*s*s*sumdy3w + 3*(s*c*c*sumdx2dyw + c*s*s*sumdxdy2w);
+    const Double_t m3t = c*c*c*sumdy3w - s*s*s*sumdx3w + 3*(s*s*c*sumdx2dyw - s*c*c*sumdxdy2w);
+
+    fM3Long  = MMath::Sqrt3(m3l/sumw);      // [mm]
+    fM3Trans = MMath::Sqrt3(m3t/sumw);      // [mm]
 
     //