Changeset 9855 for trunk/Mars
- Timestamp:
- 08/13/10 11:34:45 (14 years ago)
- Location:
- trunk/Mars
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Mars/Changelog
r9854 r9855 24 24 - added new histograms to display new variables 25 25 - increased class version number 26 27 * mimage/MHillasExt.cc: 28 - fixed claculation of the weighted time spreads 29 - fixed a bug in the calculation of the third moments! 26 30 27 31 -
trunk/Mars/mimage/MHillasExt.cc
r9850 r9855 151 151 UInt_t cnt = 0; 152 152 153 Double_t sumx = 0; 154 Double_t sumy = 0; 155 Double_t sumt = 0; 156 Double_t sumw = 0; 157 Double_t sumxy = 0; 158 Double_t sumxt = 0; 159 Double_t sumyt = 0; 160 Double_t sumxw = 0; 161 Double_t sumyw = 0; 162 Double_t sumtw = 0; 163 Double_t sumx2 = 0; 164 Double_t sumy2 = 0; 165 Double_t sumw2 = 0; 166 Double_t sumt2 = 0; 167 Double_t sumx3 = 0; 168 Double_t sumy3 = 0; 169 Double_t sumxw2 = 0; 170 Double_t sumyw2 = 0; 171 Double_t sumtx2 = 0; 172 Double_t sumty2 = 0; 173 Double_t sumtw2 = 0; 174 Double_t sumx2y = 0; 175 Double_t sumxy2 = 0; 176 Double_t sumtxw2 = 0; 177 Double_t sumtyw2 = 0; 178 Double_t sumt2w2 = 0; 153 Double_t sumx = 0; 154 Double_t sumy = 0; 155 Double_t sumt = 0; 156 Double_t sumt2 = 0; 157 Double_t sumx2 = 0; 158 Double_t sumy2 = 0; 159 Double_t sumxy = 0; 160 Double_t sumxt = 0; 161 Double_t sumyt = 0; 162 Double_t sumw = 0; 163 Double_t sumw2 = 0; 164 Double_t sumxw = 0; 165 Double_t sumyw = 0; 166 Double_t sumtw = 0; 167 Double_t sumtx2 = 0; 168 Double_t sumty2 = 0; 169 Double_t sumtw2 = 0; 170 Double_t sumxw2 = 0; 171 Double_t sumyw2 = 0; 172 Double_t sumtxw2 = 0; 173 Double_t sumtyw2 = 0; 174 Double_t sumt2w2 = 0; 175 Double_t sumx2w2 = 0; 176 Double_t sumy2w2 = 0; 177 Double_t sumxyw2 = 0; 178 179 Double_t sumdx3w = 0; 180 Double_t sumdy3w = 0; 181 182 Double_t sumdx2dyw = 0; 183 Double_t sumdxdy2w = 0; 179 184 180 185 const UInt_t npix = evt.GetNumPixels(); … … 203 208 sumt2 += t*t; 204 209 210 sumx2 += x*x; 211 sumy2 += y*y; 212 sumxy += x*y; 213 sumxt += x*t; 214 sumyt += y*t; 215 205 216 sumw += nphot; 206 217 sumw2 += nphot*nphot; … … 210 221 sumtw += t*nphot; 211 222 212 sumtx2 += x*nphot*nphot;213 sumty2 += y*nphot*nphot;223 sumtx2 += t*x*x; 224 sumty2 += t*y*y; 214 225 sumtw2 += t*nphot*nphot; 226 227 sumxw2 += x*nphot*nphot; 228 sumyw2 += y*nphot*nphot; 215 229 216 230 sumtxw2 += t*x*nphot*nphot; 217 231 sumtyw2 += t*y*nphot*nphot; 218 232 sumt2w2 += t*t*nphot*nphot; 219 220 sumx2 += x*x; 221 sumy2 += y*y; 222 sumxy += x*y; 223 sumxt += x*t; 224 sumyt += y*t; 233 sumx2w2 += x*x*nphot*nphot; 234 sumy2w2 += y*y*nphot*nphot; 235 sumxyw2 += x*y*nphot*nphot; 225 236 226 237 // --- 3rd moment --- … … 228 239 const Double_t dy = y - hil.GetMeanY(); // [mm] 229 240 230 sum x3 += nphot * dx*dx*dx;231 sum y3 += nphot * dy*dy*dy;232 sum x2y += nphot * dx*dx*dy;233 sum x2y += nphot * dx*dy*dy;241 sumdx3w += dx*dx*dx*nphot; 242 sumdy3w += dy*dy*dy*nphot; 243 sumdx2dyw += dx*dx*dy*nphot; 244 sumdxdy2w += dx*dy*dy*nphot; 234 245 235 246 cnt++; … … 253 264 254 265 // 255 // Time spread256 //257 fTimeSpread = TMath::Sqrt(sumt2*cnt - sumt*sumt)/cnt;258 fTimeSpreadWeighted = TMath::Sqrt(sumt2w2 + sumw2*sumtw*sumtw/sumw/sumw - 2*sumtw2*sumtw/sumw)/sumw;259 260 //261 266 // Time slope 262 267 // … … 277 282 278 283 // 279 // Calculate time spread around longitudinal time slope 280 // 281 const Double_t sumdt = sumt - fSlopeLong*(c*sumx + s*sumy); 282 const Double_t sumdtw = sumtw - fSlopeLong*(c*sumxw + s*sumyw); 283 const Double_t sumdtw2 = sumtw2 - fSlopeLong*(c*sumxw2 + s*sumyw2); 284 const Double_t sumdt2 = sumt2 - 2*fSlopeLong*(c*sumxt + s*sumyt) + fSlopeLong*fSlopeLong * (c*c*sumx2 + s*s*sumy2 + 2*c*s*sumxy); 285 const Double_t sumdt2w2 = sumt2w2 - 2*fSlopeLong*(c*sumtxw2 + s*sumtyw2) + fSlopeLong*fSlopeLong * (c*c*sumx2 + s*s*sumy2 + 2*c*s*sumxy) * sumw2; 286 287 fSlopeSpread = TMath::Sqrt(sumdt2*cnt - sumdt*sumdt)/cnt; 288 fSlopeSpreadWeighted = TMath::Sqrt(sumdt2w2 + sumw2*sumdtw*sumdtw/sumw/sumw - 2*sumdtw2*sumdtw/sumw)/sumw; 284 // Time spread 285 // 286 const Double_t m = fSlopeLong; 287 288 const Double_t z = c*sumx + s*sumy; 289 const Double_t zw = c*sumxw + s*sumyw; 290 const Double_t zw2 = c*sumxw2 + s*sumyw2; 291 const Double_t zt = c*sumxt + s*sumyt; 292 const Double_t ztw2 = c*sumtxw2 + s*sumtyw2; 293 294 const Double_t sumdt = sumt - m*z; 295 const Double_t sumdtw = sumtw - m*zw; 296 const Double_t sumdtw2 = sumtw2 - m*zw2; 297 298 const Double_t x2y2 = c*c*sumx2 + s*s*sumy2 + 2*c*s*sumxy; 299 const Double_t x2y2w2 = c*c*sumx2w2 + s*s*sumy2w2 + 2*c*s*sumxyw2; 300 301 const Double_t sumdt2 = sumt2 + m*m*x2y2 - 2*m*zt; 302 const Double_t sumdt2w2 = sumt2w2 + m*m*x2y2w2 - 2*m*ztw2; 303 304 const Double_t meant = sumt /cnt; 305 const Double_t meandt = sumdt/cnt; // meant - m*z/cnt 306 307 const Double_t meantw = sumtw /sumw; 308 const Double_t meandtw = sumdtw/sumw; 309 310 fTimeSpread = TMath::Sqrt(sumt2 /cnt - meant *meant); 311 fSlopeSpread = TMath::Sqrt(sumdt2/cnt - meandt*meandt); 312 313 fTimeSpreadWeighted = TMath::Sqrt(sumt2w2 + sumw2*meantw *meantw - 2*meantw *sumtw2) /sumw; 314 fSlopeSpreadWeighted = TMath::Sqrt(sumdt2w2 + sumw2*meandtw*meandtw - 2*meandtw*sumdtw2)/sumw; 289 315 290 316 // 291 317 // Third moments along axes get normalized 292 318 // 293 const Double_t m3l = c*c*c*sum x3 + 3*(s*c*c*sumx2y + c*s*s*sumxy2) + s*s*s*sumy3;294 const Double_t m3t = c*c*c*sum y3 + 3*(s*s*c*sumx2y - s*c*c*sumxy2) - s*s*s*sumx3;295 296 fM3Long = MMath::Sqrt3(m3l/ hil.GetSize()); // [mm]297 fM3Trans = MMath::Sqrt3(m3t/ hil.GetSize()); // [mm]319 const Double_t m3l = c*c*c*sumdx3w + s*s*s*sumdy3w + 3*(s*c*c*sumdx2dyw + c*s*s*sumdxdy2w); 320 const Double_t m3t = c*c*c*sumdy3w - s*s*s*sumdx3w + 3*(s*s*c*sumdx2dyw - s*c*c*sumdxdy2w); 321 322 fM3Long = MMath::Sqrt3(m3l/sumw); // [mm] 323 fM3Trans = MMath::Sqrt3(m3t/sumw); // [mm] 298 324 299 325 //
Note:
See TracChangeset
for help on using the changeset viewer.