Changeset 9860 for trunk/Mars/mimage


Ignore:
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:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Mars/mimage/MHillasExt.cc

    r9855 r9860  
    154154    Double_t sumy      = 0;
    155155    Double_t sumt      = 0;
    156     Double_t sumt2     = 0;
     156    Double_t sumw      = 0;
     157
     158    Double_t sumxt     = 0;
     159    Double_t sumyt     = 0;
     160
    157161    Double_t sumx2     = 0;
    158162    Double_t sumy2     = 0;
    159163    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 sumt2     = 0;
     165
    164166    Double_t sumxw     = 0;
    165167    Double_t sumyw     = 0;
    166168    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;
     169
     170    Double_t sumx2w    = 0;
     171    Double_t sumy2w    = 0;
     172    Double_t sumt2w    = 0;
     173    Double_t sumxyw    = 0;
     174    Double_t sumxtw    = 0;
     175    Double_t sumytw    = 0;
    178176
    179177    Double_t sumdx3w   = 0;
     
    204202        sumx    += x;
    205203        sumy    += y;
    206 
    207204        sumt    += t;
    208         sumt2   += t*t;
     205        sumw    += nphot;
     206
     207        sumxt   += x*t;
     208        sumyt   += y*t;
    209209
    210210        sumx2   += x*x;
    211211        sumy2   += y*y;
    212212        sumxy   += x*y;
    213         sumxt   += x*t;
    214         sumyt   += y*t;
    215 
    216         sumw    += nphot;
    217         sumw2   += nphot*nphot;
     213        sumt2   += t*t;
    218214
    219215        sumxw   += x*nphot;
     
    221217        sumtw   += t*nphot;
    222218
    223         sumtx2  += t*x*x;
    224         sumty2  += t*y*y;
    225         sumtw2  += t*nphot*nphot;
    226 
    227         sumxw2  += x*nphot*nphot;
    228         sumyw2  += y*nphot*nphot;
    229 
    230         sumtxw2 += t*x*nphot*nphot;
    231         sumtyw2 += t*y*nphot*nphot;
    232         sumt2w2 += t*t*nphot*nphot;
    233         sumx2w2 += x*x*nphot*nphot;
    234         sumy2w2 += y*y*nphot*nphot;
    235         sumxyw2 += x*y*nphot*nphot;
     219        sumx2w  += x*x*nphot;
     220        sumy2w  += y*y*nphot;
     221        sumt2w  += t*t*nphot;
     222
     223        sumxyw  += x*y*nphot;
     224
     225        sumxtw  += x*t*nphot;
     226        sumytw  += y*t*nphot;
    236227
    237228        // --- 3rd moment ---
     
    286277    const Double_t m = fSlopeLong;
    287278
    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 
     279    const Double_t z        = dx; //c*sumx   + s*sumy;
     280    const Double_t zw       =       c*sumxw  + s*sumyw;
     281
     282    const Double_t zt       = dxt; //c*sumxt  + s*sumyt;
     283    const Double_t ztw      =        c*sumxtw + s*sumytw;
     284
     285    const Double_t x2y2     = dx2; //c*c*sumx2  + s*s*sumy2  + 2*c*s*sumxy;
     286    const Double_t x2y2w    =        c*c*sumx2w + s*s*sumy2w + 2*c*s*sumxyw;
     287
     288    const Double_t sumdt    = sumt  - m*z;
     289    const Double_t sumdtw   = sumtw - m*zw;
     290
     291    const Double_t sumdt2   = sumt2  + m*m*x2y2  - 2*m*zt;
     292    const Double_t sumdt2w  = sumt2w + m*m*x2y2w - 2*m*ztw;
     293
     294    const Double_t meant    = sumt  /cnt;
     295    const Double_t meandt   = sumdt /cnt;
    307296    const Double_t meantw   = sumtw /sumw;
    308297    const Double_t meandtw  = sumdtw/sumw;
    309298
    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;
     299    fTimeSpread          = TMath::Sqrt(sumt2 /cnt  - meant *meant);
     300    fTimeSpreadWeighted  = TMath::Sqrt(sumt2w/sumw - meantw*meantw);
     301
     302    fSlopeSpread         = TMath::Sqrt(sumdt2 /cnt  - meandt *meandt);
     303    fSlopeSpreadWeighted = TMath::Sqrt(sumdt2w/sumw - meandtw*meandtw);
    315304
    316305    //
Note: See TracChangeset for help on using the changeset viewer.