Ignore:
Timestamp:
06/29/07 12:47:48 (17 years ago)
Author:
tbretz
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mimage/MHillasExt.cc

    r7784 r8619  
    2020!   Author(s): Wolfgang Wittek 06/2002 <mailto:wittek@mppmu.mpg.de>
    2121!
    22 !   Copyright: MAGIC Software Development, 2000-2004
     22!   Copyright: MAGIC Software Development, 2000-2007
    2323!
    2424!
    … …  
    5151// fMaxDist  added. Distance between center and most distant used pixel
    5252//
    53 //
    54 // WARNING: Before you can use fAsym, fM3Long and fM3Trans you must
     53// Version 4:
     54// ----------
     55// fMaxDist      removed.
     56// fSlopeLong    added
     57// fSlopeTrans   added
     58//
     59//
     60// WARNING: Before you can use fAsym, fSlope* and fM3* you must
    5561//          multiply by the sign of MHillasSrc::fCosDeltaAlpha
    5662//
    … …  
    6672#include <TMarker.h>
    6773#include <TVector2.h>
     74#include <TMinuit.h>
    6875#include <TVirtualPad.h>
    6976
    … …  
    97104// -------------------------------------------------------------------------
    98105//
    99 // Reset values to:
    100 //    fAsym    =  0;
    101 //    fM3Long  =  0;
    102 //    fM3Trans =  0;
    103 //    fMaxDist =  0;
     106// Reset all values to 0
    104107//
    105108void MHillasExt::Reset()
    106109{
    107     fAsym    =  0;
    108     fM3Long  =  0;
    109     fM3Trans =  0;
    110 
    111     fMaxDist =  0;
     110    fAsym       =  0;
     111    fM3Long     =  0;
     112    fM3Trans    =  0;
     113    fSlopeLong  =  0;
     114    fSlopeTrans =  0;
    112115}
    113116
    … …  
    128131    //  the complex matrix multiplication and sum is evaluated correctly.
    129132    //
    130     Double_t m3x = 0;
    131     Double_t m3y = 0;
    132 
    133133    Int_t maxpixid  = 0;
    134134    Float_t maxpix  = 0;
    135     Float_t maxdist = 0;
    136 
    137     //    MSignalPix *pix = 0;
    138     //    TIter Next(evt);
    139     //    while ((pix=(MSignalPix*)Next()))
     135
     136    // Variables to caluclate time slope
     137    // Formula: http://mo.mathematik.uni-stuttgart.de/inhalt/erlaeuterung/erlaeuterung300/
     138    UInt_t cnt = 0;
     139
     140    Double_t sumx   = 0;
     141    Double_t sumy   = 0;
     142    Double_t sumt   = 0;
     143    Double_t sumxy  = 0;
     144    Double_t sumxt  = 0;
     145    Double_t sumyt  = 0;
     146    Double_t sumx2  = 0;
     147    Double_t sumy2  = 0;
     148    Double_t sumx3  = 0;
     149    Double_t sumy3  = 0;
     150    Double_t sumx2y = 0;
     151    Double_t sumxy2 = 0;
    140152
    141153    const UInt_t npix = evt.GetNumPixels();
    … …  
    149161            continue;
    150162
    151         //const Int_t pixid = pix->GetPixId();
    152 
    153         const MGeomPix &gpix = geom[i/*pixid*/];
    154         const Double_t dx = gpix.GetX() - hil.GetMeanX();      // [mm]
    155         const Double_t dy = gpix.GetY() - hil.GetMeanY();      // [mm]
    156 
    157         Double_t nphot = pix.GetNumPhotons();                  // [1]
    158 
    159         const Double_t dzx =  hil.GetCosDelta()*dx + hil.GetSinDelta()*dy; // [mm]
    160         const Double_t dzy = -hil.GetSinDelta()*dx + hil.GetCosDelta()*dy; // [mm]
    161 
    162         const Double_t dist = dx*dx+dy*dy;
    163         if (dist>TMath::Abs(maxdist))
    164             maxdist = dzx<0 ? -dist : dist;                    // [mm^2]
    165 
    166         m3x += nphot * dzx*dzx*dzx;                            // [mm^3]
    167         m3y += nphot * dzy*dzy*dzy;                            // [mm^3]
     163        const MGeomPix &gpix = geom[i];
     164
     165        const Double_t x = gpix.GetX();
     166        const Double_t y = gpix.GetY();
     167        const Double_t t = pix.GetArrivalTime();
     168
     169        // --- time slope ----
     170        sumx  += x;
     171        sumy  += y;
     172        sumt  += t;
     173
     174        sumx2 += x*x;
     175        sumy2 += y*y;
     176        sumxy += x*y;
     177        sumxt += x*t;
     178        sumyt += y*t;
     179
     180        // --- 3rd moment ---
     181        const Double_t dx = x - hil.GetMeanX();      // [mm]
     182        const Double_t dy = y - hil.GetMeanY();      // [mm]
     183
     184        Double_t nphot = pix.GetNumPhotons();        // [1]
     185
     186        sumx3  += nphot * dx*dx*dx;
     187        sumy3  += nphot * dy*dy*dy;
     188        sumx2y += nphot * dx*dx*dy;
     189        sumx2y += nphot * dx*dy*dy;
     190
     191        cnt++;
    168192
    169193        //
    … …  
    171195        // must take pixel size into account
    172196        //
    173         nphot *= geom.GetPixRatio(i/*pixid*/);
    174 
     197        nphot *= geom.GetPixRatio(i);
     198
     199        // --- max pixel ---
    175200        if (nphot>maxpix)
    176201        {
    177             maxpix   = nphot;                                  // [1]
    178             maxpixid = i;//pixid;
     202            maxpix   = nphot;                        // [1]
     203            maxpixid = i;
    179204        }
    180205    }
    181206
     207    const Double_t c = hil.GetCosDelta();
     208    const Double_t s = hil.GetSinDelta();
     209
     210    //
     211    // Time slope
     212    //
     213    const Double_t dxt  =  c*sumxt + s*sumyt;
     214    const Double_t dyt  = -s*sumxt + c*sumyt;
     215
     216    const Double_t dx   =  c*sumx  + s*sumy;
     217    const Double_t dy   = -s*sumx  + c*sumy;
     218
     219    const Double_t dx2  =  c*c*sumx2 + 2*c*s*sumxy + s*s*sumy2;
     220    const Double_t dy2  =  s*s*sumx2 - 2*c*s*sumxy + c*c*sumy2;
     221
     222    const Double_t detx =  cnt*dx2 - dx*dx;
     223    const Double_t dety =  cnt*dy2 - dy*dy;
     224
     225    fSlopeLong  = detx==0 ? 0 : (cnt*dxt - sumt*dx)/detx;
     226    fSlopeTrans = dety==0 ? 0 : (cnt*dyt - sumt*dy)/dety;
     227
     228    //
     229    // Third moments along axes get normalized
     230    //
     231    const Double_t m3l = c*c*c*sumx3 + 3*(s*c*c*sumx2y + c*s*s*sumxy2) + s*s*s*sumy3;
     232    const Double_t m3t = c*c*c*sumy3 + 3*(s*s*c*sumx2y - s*c*c*sumxy2) - s*s*s*sumx3;
     233
     234    fM3Long  = MMath::Sqrt3(m3l/hil.GetSize());      // [mm]
     235    fM3Trans = MMath::Sqrt3(m3t/hil.GetSize());      // [mm]
     236
     237    //
     238    // Asymmetry
     239    //
    182240    const MGeomPix &maxp = geom[maxpixid];
    183 
    184     //
    185     // Asymmetry
    186     //
    187     fAsym = (hil.GetMeanX()-maxp.GetX())*hil.GetCosDelta() +
    188             (hil.GetMeanY()-maxp.GetY())*hil.GetSinDelta();            // [mm]
    189 
    190     //
    191     // Third moments along axes get normalized
    192     //
    193     m3x /= hil.GetSize();
    194     m3y /= hil.GetSize();
    195 
    196     fM3Long  = m3x<0 ? -pow(-m3x, 1./3) : pow(m3x, 1./3);      // [mm]
    197     fM3Trans = m3y<0 ? -pow(-m3y, 1./3) : pow(m3y, 1./3);      // [mm]
    198 
    199     //
    200     // Distance between max signal and COG
    201     //
    202     const Double_t md = TMath::Sqrt(TMath::Abs(maxdist));
    203     fMaxDist = maxdist<0 ? -md : md;                           // [mm]
     241    fAsym = (hil.GetMeanX()-maxp.GetX())*c + (hil.GetMeanY()-maxp.GetY())*s;            // [mm]
    204242
    205243    SetReadyToSave();
    … …  
    215253void MHillasExt::Set(const TArrayF &arr)
    216254{
    217     if (arr.GetSize() != 4)
     255    if (arr.GetSize() != 5)
    218256        return;
    219257
    220     fAsym    = arr.At(0);
    221     fM3Long  = arr.At(1);
    222     fM3Trans = arr.At(2);
    223     fMaxDist = arr.At(3);
     258    fAsym       = arr.At(0);
     259    fM3Long     = arr.At(1);
     260    fM3Trans    = arr.At(2);
     261    fSlopeLong  = arr.At(3);
     262    fSlopeTrans = arr.At(4);
    224263}
    225264
    … …  
    228267// Print contents of MHillasExt to *fLog.
    229268//
    230 void MHillasExt::Print(Option_t *) const
    231 {
     269void MHillasExt::Print(Option_t *o) const
     270{
     271    const Bool_t showtrans = TString(o).Contains("trans");
     272
    232273    *fLog << all;
    233274    *fLog << GetDescriptor() << endl;
    234275    *fLog << " - Asymmetry      [mm]  = " << fAsym    << endl;
    235276    *fLog << " - 3.Moment Long  [mm]  = " << fM3Long  << endl;
    236     *fLog << " - 3.Moment Trans [mm]  = " << fM3Trans << endl;
    237     *fLog << " - Max.Dist       [mm]  = " << fMaxDist << endl;
     277    if (showtrans)
     278        *fLog << " - 3.Moment Trans [mm]  = " << fM3Trans << endl;
     279    *fLog << " - Slope Long     [mm]  = " << fSlopeLong  << endl;
     280    if (showtrans)
     281        *fLog << " - Slope Trans    [mm]  = " << fSlopeTrans << endl;
    238282}
    239283
    … …  
    249293    *fLog << " - Asymmetry      [deg] = " << fAsym   *geom.GetConvMm2Deg() << endl;
    250294    *fLog << " - 3.Moment Long  [deg] = " << fM3Long *geom.GetConvMm2Deg() << endl;
    251     *fLog << " - 3.Moment Trans [deg] = " << fM3Trans*geom.GetConvMm2Deg() << endl;
    252     *fLog << " - Max.Dist       [deg] = " << fMaxDist*geom.GetConvMm2Deg() << endl;
     295    *fLog << " - 3.Moment Trans [mm]  = " << fM3Trans*geom.GetConvMm2Deg() << endl;
     296    *fLog << " - Slope Long     [mm]  = " << fSlopeLong*geom.GetConvMm2Deg()  << endl;
     297    *fLog << " - Slope Trans    [mm]  = " << fSlopeTrans*geom.GetConvMm2Deg() << endl;
    253298}
    254299
    … …  
    272317
    273318    TMarker m;
    274     m.SetMarkerColor(15);
     319    m.SetMarkerColor(kRed);
    275320    m.SetMarkerStyle(kFullDotLarge);
    276321    m.PaintMarker(v.X(), v.Y());
Note: See TracChangeset for help on using the changeset viewer.