Changeset 9855 for trunk/Mars


Ignore:
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!
Location:
trunk/Mars
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/Mars/Changelog

    r9854 r9855  
    2424     - added new histograms to display new variables
    2525     - 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!
    2630
    2731
  • trunk/Mars/mimage/MHillasExt.cc

    r9850 r9855  
    151151    UInt_t cnt = 0;
    152152
    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;
    179184
    180185    const UInt_t npix = evt.GetNumPixels();
     
    203208        sumt2   += t*t;
    204209
     210        sumx2   += x*x;
     211        sumy2   += y*y;
     212        sumxy   += x*y;
     213        sumxt   += x*t;
     214        sumyt   += y*t;
     215
    205216        sumw    += nphot;
    206217        sumw2   += nphot*nphot;
     
    210221        sumtw   += t*nphot;
    211222
    212         sumtx2  += x*nphot*nphot;
    213         sumty2  += y*nphot*nphot;
     223        sumtx2  += t*x*x;
     224        sumty2  += t*y*y;
    214225        sumtw2  += t*nphot*nphot;
     226
     227        sumxw2  += x*nphot*nphot;
     228        sumyw2  += y*nphot*nphot;
    215229
    216230        sumtxw2 += t*x*nphot*nphot;
    217231        sumtyw2 += t*y*nphot*nphot;
    218232        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;
    225236
    226237        // --- 3rd moment ---
     
    228239        const Double_t dy = y - hil.GetMeanY();      // [mm]
    229240
    230         sumx3  += nphot * dx*dx*dx;
    231         sumy3  += nphot * dy*dy*dy;
    232         sumx2y += nphot * dx*dx*dy;
    233         sumx2y += 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;
    234245
    235246        cnt++;
     
    253264
    254265    //
    255     // Time spread
    256     //
    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     //
    261266    // Time slope
    262267    //
     
    277282
    278283    //
    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;
    289315
    290316    //
    291317    // Third moments along axes get normalized
    292318    //
    293     const Double_t m3l = c*c*c*sumx3 + 3*(s*c*c*sumx2y + c*s*s*sumxy2) + s*s*s*sumy3;
    294     const Double_t m3t = c*c*c*sumy3 + 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]
    298324
    299325    //
Note: See TracChangeset for help on using the changeset viewer.