Index: trunk/Mars/Changelog =================================================================== --- trunk/Mars/Changelog (revision 9854) +++ trunk/Mars/Changelog (revision 9855) @@ -24,4 +24,8 @@ - added new histograms to display new variables - increased class version number + + * mimage/MHillasExt.cc: + - fixed claculation of the weighted time spreads + - fixed a bug in the calculation of the third moments! Index: trunk/Mars/mimage/MHillasExt.cc =================================================================== --- trunk/Mars/mimage/MHillasExt.cc (revision 9854) +++ trunk/Mars/mimage/MHillasExt.cc (revision 9855) @@ -151,30 +151,35 @@ 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,4 +208,10 @@ sumt2 += t*t; + sumx2 += x*x; + sumy2 += y*y; + sumxy += x*y; + sumxt += x*t; + sumyt += y*t; + sumw += nphot; sumw2 += nphot*nphot; @@ -210,17 +221,17 @@ 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,8 +239,8 @@ 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,10 +264,4 @@ // - // 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,23 +282,44 @@ // - // 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] //