Changeset 8619 for trunk/MagicSoft


Ignore:
Timestamp:
06/29/07 12:47:48 (17 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r8618 r8619  
    2020
    2121
     22 2007/06/29 Thomas Bretz
     23
     24   * mimage/MHillasExt.[h,cc]:
     25     - added new timing parameters fSlopeTrans and fSlopeLong
     26     - removed never used parameter fMaxDist
     27     - increased accordingly the class version number by 1
     28     - replaced the old calculation of the third moments by a
     29       numerically more accurate one, by avoiding to calculate
     30       to many differences too often.
     31
     32   * mfilter/MFMagicCuts.[h,cc]:
     33     - added the usage of the new SlopeLong parameter for ghostbusting
     34
     35
     36
    2237 2007/06/28 Thomas Bretz
    2338
    … …  
    4863   * mjobs/MJCalibrateSignal.cc:
    4964     - changed axis label of PulsePos plot (now in nanosec)
     65
     66   * mpointing/MSrcPosCalc.[h,cc]:
     67     - allow to set a tasklist as callback to now which n-th
     68     pass of the same task list it is
     69
     70   * mbase/MTaskList.[h,cc]:
     71     - added some code to allow the execution of one task list more
     72       than once. This is for example necessary to process three
     73       different off-source regions.
     74
     75   * mjobs/MJCut.[h,cc]:
     76     - use the new feature in MTaskList to setup a tasklist
     77       processing the off-source calculation tasklist more than once
     78     - added a new data meber fNumOffSourcePos
     79     - added a new resource option NumOffSourcePositions
     80     - added a new CutQ before Cut0 which takes place before all
     81       source posisiton dependant stuff
     82
     83   * ganymed_onoff.rc:
     84     - renamed Cut0 to CutQ
    5085
    5186

  • trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc

    r8613 r8619  
    1818!   Author(s): Thomas Bretz, 03/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
    1919!
    20 !   Copyright: MAGIC Software Development, 2000-2005
     20!   Copyright: MAGIC Software Development, 2000-2007
    2121!
    2222!
    … …  
    3838//      dist         = MHillasSrc.fDist*fMm2Deg
    3939//      m3long       = MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*fMm2Deg
     40//      slope        = MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)*fMm2Deg
    4041//
    4142//      antialpha    = MHillasSrcAnti.fAlpha
    4243//      antidist     = MHillasSrcAnti.fDist*fMm2Deg
    4344//      antim3long   = MHillasExt.fM3Long*sign(MHillasSrcAnti.fCosDeltaAlpha)*fMm2Deg
     45//      antislope    = MHillasExt.fSlopeLong*sign(MHillasSrcAnti.fCosDeltaAlpha)*fMm2Deg
    4446//
    4547//      had          = Hadronness.fVal
    … …  
    4850//  ------------------
    4951//
    50 //    c[0], c[5], c[6], c[7]:
     52//    c[0],  c[5], c[6], c[7],  c[8], c[9]:
    5153//          xi          = c[0]+c[6]*pow(leakage1, c[7]);
    5254//          p           = xi*(1-width/length);
    53 //          disp        = TMath::Sign(disp, m3long-c[5])
    54 //          antidisp    = TMath::Sign(disp, antim3long-c[5])
     55//          sign1       = m3long-c[5]
     56//          sign2       = (dist-c[9])*c[8]-slope
     57//          disp        = sign1<0 ||sign2<0 ? -disp : disp
     58//          antisign1   = antim3long-c[5]
     59//          antisign2   = (antidist-c[9])*c[8]-antislope
     60//          antidisp    = antisign1<0 || antisign2<0 ? -antidisp : antidisp
    5561//          thetasq     = disp^2 + dist^2 - 2*disp*dist*alpha
    5662//          antithetasq = antidisp^2 + antidist^2 - 2*antidisp*antidist*antialpha
    … …  
    6571//
    6672//    HadronnessCut:
    67 //    c[8], c[9]:
    68 //          had       <= c[8]
    69 //          10^lgsize >= c[9]
     73//    c[10], c[11]:
     74//          had       <= c[10]
     75//          10^lgsize >= c[11]
    7076//
    7177//
    … …  
    150156    AddToBranchList("MHillasSrcAnti.fCosDeltaAlpha");
    151157    AddToBranchList("MHillasExt.fM3Long");
     158    AddToBranchList("MHillasExt.fSlopeLong");
    152159    AddToBranchList("MNewImagePar.fLeakage1");
    153160    AddToBranchList("Hadronness.fVal");
    … …  
    277284    fMatrix = mat;
    278285
     286//    fMap[kELength] = fMatrix->AddColumn("MHillas.fLength*MGeomCam.fConvMm2Deg");
     287//    fMap[kEWidth]  = fMatrix->AddColumn("MHillas.fWidth*MGeomCam.fConvMm2Deg");
     288
    279289    fMap[kEWdivL]   = fMatrix->AddColumn("MHillas.fWidth/MHillas.fLength");
    280290    fMap[kESize]    = fMatrix->AddColumn("log10(MHillas.fSize)");
    … …  
    289299    fMap[kESrcSign] = fMatrix->AddColumn("sign(MHillasSrc.fCosDeltaAlpha)");
    290300
     301    fMap[kESlope]   = fMatrix->AddColumn("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/MGeomCam.fConvMm2Deg");
     302
    291303    if (fThetaCut&kOff)
    292304    {
    … …  
    294306        fMap[kEDistAnti]   = fMatrix->AddColumn("MHillasSrcAnti.fDist*MGeomCam.fConvMm2Deg");
    295307        fMap[kEM3LongAnti] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrcAnti.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");
     308        fMap[kESlopeAnti]  = fMatrix->AddColumn("MHillasExt.fSlopeLong*sign(MHillasSrcAnti.fCosDeltaAlpha)/MGeomCam.fConvMm2Deg");
    296309    }
    297310
    … …  
    351364    const Double_t dist   = fMatrix ? GetVal(kEDist)    : fHilSrc->GetDist()*fMm2Deg;
    352365    const Double_t alpha  = fMatrix ? GetVal(kEAlpha)   : fHilSrc->GetAlpha();
     366    const Double_t sign   = fMatrix ? GetVal(kESrcSign) : MMath::Sgn(fHilSrc->GetCosDeltaAlpha());
    353367    const Double_t m3long = fMatrix ? GetVal(kEM3Long)  : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha());
    354     const Double_t sign   = fMatrix ? GetVal(kESrcSign) : MMath::Sgn(fHilSrc->GetCosDeltaAlpha());
     368    const Double_t slope  = fMatrix ? GetVal(kESlope)   : fHilExt->GetSlopeLong()/TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha());
     369
     370        // Do Ghostbusting with slope and m3l
     371    const Double_t sign1 = (dist-c[9])*c[8]-slope;
     372    const Double_t sign2 = m3long-c[5];
     373    const Double_t p = sign1<0 || sign2<0 ? -disp : disp;
    355374
    356375    // Align disp along source direction depending on third moment
    357     const Double_t p = TMath::Sign(disp, m3long-c[5]);
     376    //const Double_t p = TMath::Sign(disp, m3long-c[5]);
    358377
    359378    // Align sign of disp along shower axis like m3long
    … …  
    384403    {
    385404        const Double_t had = fMatrix ? GetVal(kEHadronness) : fHadronness->GetVal();
    386         if (had>c[8])
     405        if (had>c[10])
    387406            return kTRUE;
    388407
    389         if (TMath::Power(10, lgsize)<c[9])
     408        if (TMath::Power(10, lgsize)<c[11])
    390409            return kTRUE;
    391410    }
    … …  
    399418        const Double_t alphaanti   = fMatrix ? GetVal(kEAlphaAnti)  : fHilAnti->GetAlpha();
    400419        const Double_t m3lanti     = fMatrix ? GetVal(kEM3LongAnti) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha());
    401 
    402         const Double_t panti       = TMath::Sign(disp, m3lanti-c[5]);
     420        const Double_t slopeanti   = fMatrix ? GetVal(kESlopeAnti)  : fHilExt->GetSlopeLong()/TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha());
     421
     422        // Do Ghostbusting with slope and m3l
     423        const Double_t sign3 = (distanti-c[9])*c[8]-slopeanti;
     424        const Double_t sign4 = m3lanti-c[5];
     425        const Double_t panti = sign3<0 || sign4<0 ? -disp : disp;
     426
     427        // Align disp along source direction depending on third moment
     428        //const Double_t panti       = TMath::Sign(disp, m3lanti-c[5]);
    403429        const Double_t thetasqanti = GetThetaSq(panti, wdivl, distanti, alphaanti);
    404430

  • trunk/MagicSoft/Mars/mfilter/MFMagicCuts.h

    r7979 r8619  
    4343    enum {
    4444        kESize, kEAlpha, kEAlphaAnti, kEArea, kEDist, kEDistAnti,
    45         kEM3Long, kEM3LongAnti, kEWdivL, kELeakage, kESrcSign, kEHadronness,
     45        kEM3Long, kEM3LongAnti, kEWdivL, kELeakage, kESrcSign,
     46        kESlope, kESlopeAnti, kEHadronness,
    4647        kLastElement
    4748    };

  • 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());

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

    r7176 r8619  
    2020    Float_t fM3Trans; // [mm] 3rd moment (e-weighted) along minor axis
    2121
    22     Float_t fMaxDist; // Distance between center and most distant used pixel
     22    Float_t fSlopeLong;
     23    Float_t fSlopeTrans;
    2324
    2425public:
    … …  
    2728    void Reset();
    2829
    29     Float_t GetAsym() const    { return fAsym; }
    30     Float_t GetM3Long() const  { return fM3Long; }
    31     Float_t GetM3Trans() const { return fM3Trans; }
    32 
    33     Float_t GetMaxDist() const { return fMaxDist; }
     30    Float_t GetAsym() const       { return fAsym; }
     31    Float_t GetM3Long() const     { return fM3Long; }
     32    Float_t GetM3Trans() const    { return fM3Trans; }
     33    Float_t GetSlopeLong() const  { return fSlopeLong; }
     34    Float_t GetSlopeTrans() const { return fSlopeTrans; }
    3435
    3536    Int_t Calc(const MGeomCam &geom, const MSignalCam &pix,
    … …  
    4344    void Set(const TArrayF &arr);
    4445
    45     ClassDef(MHillasExt, 3) // Storage Container for extended Hillas Parameter
     46    ClassDef(MHillasExt, 4) // Storage Container for extended Hillas Parameter
    4647};
    4748#endif
Note: See TracChangeset for help on using the changeset viewer.