Ignore:
Timestamp:
07/14/05 17:13:40 (20 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mfilter
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r7173 r7191  
    2727//   MFMagicCuts
    2828//
     29//  Predefinitions:
     30//  ---------------
     31//
     32//      m3long     = MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*fMM2Deg
     33//      antim3long = MHillasExt.fM3Long*sign(MHillasSrcAnti.fCosDeltaAlpha)*fMM2Deg
     34//
     35//  Coefficients/Cuts:
     36//  ------------------
     37//
     38//    c[0], c[5], c[6], c[7]:
     39//          xi          = c[0]+c[6]*pow(Leakage1, c[7]);
     40//          p           = xi*(1-width/length);
     41//          disp        = TMath::Sign(disp, m3long-c[5])
     42//          antidisp    = TMath::Sign(disp, antim3long-c[5])
     43//          thetasq     = disp^2 + dist^2 - 2*disp*dist*alpha
     44//          antithetasq = antidisp^2 + dist^2 - 2*antidisp*dist*alpha
     45//
     46//    c[1]:
     47//          thetasq < c[1]*c[1]
     48//
     49//    c[2], c[3], c[4]:
     50//          A = c[2]*(1 - c[4]*(lgsize-c[3])*(lgsize-c[3]))
     51//          area < A
     52//
     53//
     54//  Options:
     55//  --------
     56//
     57//   HadronnessCut:
     58//    - none/nocut: No area cut
     59//    - area:       Area cut
     60//    - all:        same as area
     61//
     62//   ThetaCut:
     63//    - none/nocut: no theta cut
     64//    - on:         cut in theta only
     65//    - off:        anti-cut in anti-theta only
     66//    - wobble:     same as on|off (both)
     67//
     68//
     69//  Input Container:
     70//  ------
     71//   MGeomCam
     72//   MHillas
     73//   MHillasExt
     74//   MNewImagePar
     75//   MHillasSrc
     76//   [MHillasSrcAnti [MHillasSrc]]
     77//
     78//  Output:
     79//  -------
     80//   ThetaSquared [MParameterD]    // stores thetasq
     81//   Disp [MParameterD]            // stores -p*sign(MHillasSrc.fCosDeltaAlpha)
     82//   ThetaSquaredCut [MParameterD] // stores c[1]*c[1]
     83//
    2984/////////////////////////////////////////////////////////////////////////////
    3085#include "MFMagicCuts.h"
    … …  
    3691#include "MLog.h"
    3792#include "MLogManip.h"
     93
     94#include "MMath.h"
    3895
    3996#include "MParList.h"
    … …  
    4299#include "MHillasSrc.h"
    43100#include "MHillasExt.h"
     101#include "MNewImagePar.h"
    44102#include "MParameters.h"
    45 #include "MPointingPos.h"
    46103
    47104ClassImp(MFMagicCuts);
    … …  
    54111//
    55112MFMagicCuts::MFMagicCuts(const char *name, const char *title)
    56     : fHil(0), fHilSrc(0), fHilAnti(0), fHilExt(0), fThetaSq(0),
    57     fDisp(0), fPointing(0), fMatrix(0), fVariables(30), fThetaCut(kNone),
     113    : fHil(0), fHilSrc(0), fHilAnti(0), fHilExt(0), fNewImgPar(0),
     114    fThetaSq(0), fDisp(0), fMatrix(0), fVariables(30), fThetaCut(kNone),
    58115    fHadronnessCut(kAll)
    59116{
    … …  
    68125    AddToBranchList("MHillasSrc.fAlpha");
    69126    AddToBranchList("MHillasSrcAnti.fAlpha");
    70     AddToBranchList("MPointingPos.fZd");
    71127    AddToBranchList("MHillasExt.fMaxDist");
    72128    AddToBranchList("MHillasExt.fM3Long");
     129    AddToBranchList("MNewImagePar.fLeakage1");
    73130
    74131    fVariables[0] =  1.42547;      // Xi
    … …  
    110167        return kFALSE;
    111168    thetacut->SetVal(GetThetaSqCut());
    112     //thetacut->SetReadyToSave();
    113 
    114     MParameterD *m3lcut = (MParameterD*)pList->FindCreateObj("MParameterD", "M3lCut");
    115     if (!m3lcut)
    116         return kFALSE;
    117     m3lcut->SetVal(fVariables[5]);
    118 
    119     MParameterD *dispxi = (MParameterD*)pList->FindCreateObj("MParameterD", "DispXi");
    120     if (!dispxi)
    121         return kFALSE;
    122     dispxi->SetVal(fVariables[0]);
    123 
    124     dispxi = (MParameterD*)pList->FindCreateObj("MParameterD", "DispXiTheta");
    125     if (!dispxi)
    126         return kFALSE;
    127     dispxi->SetVal(fVariables[6]);
    128169
    129170    Print();
    … …  
    148189    }
    149190
    150     fPointing = (MPointingPos*)pList->FindObject("MPointingPos");
    151     if (!fPointing)
    152     {
    153         *fLog << err << "MPointingPos not found... aborting." << endl;
     191    fNewImgPar = (MNewImagePar*)pList->FindObject("MNewImagePar");
     192    if (!fNewImgPar)
     193    {
     194        *fLog << err << "MNewImagePar not found... aborting." << endl;
    154195        return kFALSE;
    155196    }
    … …  
    200241    fMatrix = mat;
    201242
    202     fMap[kESize]   = fMatrix->AddColumn("log10(MHillas.fSize)");
    203     fMap[kEArea]   = fMatrix->AddColumn("MHillas.GetArea*MGeomCam.fConvMm2Deg*MGeomCam.fConvMm2Deg");
    204     fMap[kEM3Long] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");
    205     fMap[kEWdivL]  = fMatrix->AddColumn("MHillas.fWidth/MHillas.fLength");
    206     fMap[kEZd]     = fMatrix->AddColumn("MPointingPos.GetZdRad");
    207 
    208     fMap[kEAlpha]  = fMatrix->AddColumn("MHillasSrc.fAlpha");
    209     fMap[kEDist]   = fMatrix->AddColumn("MHillasSrc.fDist*MGeomCam.fConvMm2Deg");
    210     //fMap[kDCA]    = fMatrix->AddColumn("MHillasSrc.fDCA*MGeomCam.fConvMm2Deg");
     243    fMap[kESize]    = fMatrix->AddColumn("log10(MHillas.fSize)");
     244    fMap[kEArea]    = fMatrix->AddColumn("MHillas.GetArea*MGeomCam.fConvMm2Deg*MGeomCam.fConvMm2Deg");
     245    fMap[kEM3Trans] = fMatrix->AddColumn("abs(MHillasExt.fM3Trans)*MGeomCam.fConvMm2Deg");
     246    fMap[kEM3Long]  = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");
     247    fMap[kESrcSign] = fMatrix->AddColumn("sign(MHillasSrc.fCosDeltaAlpha)");
     248    fMap[kEWdivL]   = fMatrix->AddColumn("MHillas.fWidth/MHillas.fLength");
     249
     250    fMap[kEAlpha]   = fMatrix->AddColumn("MHillasSrc.fAlpha");
     251    fMap[kEDist]    = fMatrix->AddColumn("MHillasSrc.fDist*MGeomCam.fConvMm2Deg");
     252    fMap[kELeakage] = fMatrix->AddColumn("log10(MNewImagePar.fLeakage1+1)");
     253
    211254    if (fThetaCut&kOff)
    212255    {
    213256        fMap[kEAlphaAnti]  = fMatrix->AddColumn("MHillasSrcAnti.fAlpha");
    214257        fMap[kEDistAnti]   = fMatrix->AddColumn("MHillasSrcAnti.fDist*MGeomCam.fConvMm2Deg");
    215         fMap[kEM3LongAnti] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrcAnti.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");
    216         //fMap[kDCAAnti]   = fMatrix->AddColumn("MHillasSrcAnti.fDCA*MGeomCam.fConvMm2Deg");
    217     }
    218 
    219     //fMap[kLength] = fMatrix->AddColumn("MHillas.fLength*MGeomCam.fConvMm2Deg");
    220     //fMap[kWidth]  = fMatrix->AddColumn("MHillas.fWidth*MGeomCam.fConvMm2Deg");
     258    }
    221259}
    222260
    … …  
    229267// If par!=NULL p is stored in this MParameterD
    230268//
    231 Double_t MFMagicCuts::GetThetaSq(Double_t c, Double_t wl, Double_t d, Double_t a, MParameterD *par) const
     269Double_t MFMagicCuts::GetThetaSq(Double_t p, Double_t wl, Double_t d, Double_t a) const
    232270{
    233271    // wl = width/l [1]
    234272    // d  = dist    [deg]
    235273    // a  = alpha   [deg]
    236     const Double_t p = c*(1-wl);
    237     if (par)
    238         par->SetVal(p);
    239     return d*d + p*p - 2*d*p*cos(a*TMath::DegToRad());
     274    return d*d + p*p - 2*d*p*TMath::Cos(a*TMath::DegToRad());
    240275}
    241276
    … …  
    256291{
    257292    // Get some variables
    258     const Double_t wdivl  = fMatrix ? GetVal(kEWdivL) : fHil->GetWidth()/fHil->GetLength();
    259     const Double_t lgsize = fMatrix ? GetVal(kESize)  : TMath::Log10(fHil->GetSize());
    260     const Double_t zd     = fMatrix ? GetVal(kEZd)    : fPointing->GetZdRad();
     293    const Double_t wdivl  = fMatrix ? GetVal(kEWdivL)   : fHil->GetWidth()/fHil->GetLength();
     294    const Double_t lgsize = fMatrix ? GetVal(kESize)    : TMath::Log10(fHil->GetSize());
     295    //const Double_t zd     = fMatrix ? GetVal(kEZd)      : fPointing->GetZdRad();
     296    const Double_t leak   = fMatrix ? GetVal(kELeakage) : TMath::Log10(fNewImgPar->GetLeakage1()+1);
    261297
    262298    // For simplicity
    263299    const Double_t *c = fVariables.GetArray();
    264300
    265     // Value for Theta cut
    266     const Double_t cut = GetThetaSqCut(); //c[1]*c[1];
    267     const Double_t xi  = (c[0]+c[8]*(lgsize-c[7])*(lgsize-c[7])) - c[6]*zd*zd;
     301    // Value for Theta cut (Disp parametrization)
     302    const Double_t cut  = GetThetaSqCut();
     303    const Double_t xi   = c[0]+c[6]*pow(leak, c[7]);
     304    const Double_t disp = xi*(1-wdivl);
    268305
    269306    // Default if we return before the end
    … …  
    272309    // ------------------- Most efficient cut -----------------------
    273310    // ---------------------- Theta cut ON --------------------------
    274     const Double_t dist    = fMatrix ? GetVal(kEDist)  : fHilSrc->GetDist()*fMm2Deg;
    275     const Double_t alpha   = fMatrix ? GetVal(kEAlpha) : fHilSrc->GetAlpha();
    276 
    277     const Double_t thetasq = GetThetaSq(xi, wdivl, dist, alpha, fDisp);
    278 
     311    const Double_t dist   = fMatrix ? GetVal(kEDist)    : fHilSrc->GetDist()*fMm2Deg;
     312    const Double_t alpha  = fMatrix ? GetVal(kEAlpha)   : fHilSrc->GetAlpha();
     313    const Double_t m3long = fMatrix ? GetVal(kEM3Long)  : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha());
     314    const Double_t sign   = fMatrix ? GetVal(kESrcSign) : MMath::Sgn(fHilSrc->GetCosDeltaAlpha());
     315
     316    // Align disp along source direction depending on third moment
     317    const Double_t p = TMath::Sign(disp, m3long-c[5]);
     318
     319    // Align sign of disp along shower axis like m3long
     320    fDisp->SetVal(-p*sign);
     321
     322    // Calculate corresponding Theta^2
     323    const Double_t thetasq = GetThetaSq(p, wdivl, dist, alpha);
    279324    fThetaSq->SetVal(thetasq);
    280325
     326    // Perform theta cut
    281327    if (fThetaCut&kOn && thetasq >= cut)
    282328        return kTRUE;
    … …  
    290336        if (area>=A)
    291337            return kTRUE;
    292     }
    293     if (fHadronnessCut&kM3Long)
    294     {
    295         const Double_t m3long = fMatrix ? GetVal(kEM3Long) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha());
    296         if (m3long<=c[5])
    297             return kTRUE;
     338
     339        //const Double_t m3t  = fMatrix ? GetVal(kEM3Trans) : TMath::Abs(fHilExt->GetM3Trans())*fMm2Deg;
     340        //if (m3t>c[8])
     341        //    return kTRUE;
    298342    }
    299343
    300344    // ------------------- Least efficient cut -----------------------
    301345    // ---------------------- Theta cut OFF --------------------------
     346
    302347    if (fThetaCut&kOff)
    303348    {
    304349        const Double_t distanti    = fMatrix ? GetVal(kEDistAnti)   : fHilAnti->GetDist()*fMm2Deg;
    305350        const Double_t alphaanti   = fMatrix ? GetVal(kEAlphaAnti)  : fHilAnti->GetAlpha();
    306         const Double_t thetasqanti = GetThetaSq(xi, wdivl, distanti, alphaanti);
    307351        const Double_t m3lanti     = fMatrix ? GetVal(kEM3LongAnti) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha());
    308352
    309         if (thetasqanti<cut && (fHadronnessCut&kM3Long || m3lanti>c[5]))
     353        const Double_t panti       = TMath::Sign(disp, m3lanti-c[5]);
     354        const Double_t thetasqanti = GetThetaSq(panti, wdivl, distanti, alphaanti);
     355
     356        if (thetasqanti<cut)
    310357            return kTRUE;
     358
     359        // This cut throws away gamma-like events which would be assigned to the
     360        // anti-source. It is not necessary but it offers the possibility to
     361        // fill the MHDisp plot in one run (Having a hole instead of eg. Crab
     362        // the data can be rotated and subtracted)
    311363    }
    312364
    … …  
    355407        *fLog << "none" << endl;
    356408        break;
    357     case kArea:
    358         *fLog << "area" << endl;
    359         break;
    360     case kM3Long:
    361         *fLog << "m3long" << endl;
    362         break;
     409    //case kArea:
     410    //    *fLog << "area" << endl;
     411    //    break;
    363412    case kAll:
    364413        *fLog << "all" << endl;
    … …  
    449498        if (str==(TString)"area")
    450499            fHadronnessCut = kArea;
    451         if (str==(TString)"m3long")
    452             fHadronnessCut = kM3Long;
    453500        if (str==(TString)"all")
    454501            fHadronnessCut = kAll;

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

    r7173 r7191  
    1515class MHillasSrc;
    1616class MHillasExt;
     17class MNewImagePar;
    1718class MParameterD;
    1819class MPointingPos;
    … …  
    3334        kNoCut =BIT(0),
    3435        kArea  =BIT(1),
    35         kM3Long=BIT(2),
    36         kAll   =kArea|kM3Long
     36        kAll   =kArea
    3737    };
    3838
    … …  
    4141    // the last on in the list. It is used as counter for fMap.
    4242    enum {
    43         kESize, kEAlpha, kEAlphaAnti, kEArea, kEDist,
    44         kEM3Long, kEM3LongAnti, kEDistAnti, kEWdivL, kEZd,
     43        kESize, kEAlpha, kEAlphaAnti, kEArea, kEDist, kEM3Long,
     44        kEM3LongAnti, kEM3Trans, kEDistAnti, kEWdivL, //kEZd,
     45        kELeakage, kESrcSign, //kEDelta, //kEMeanX, kEMeanY, kEDelta,
    4546        kLastElement
    4647    };
    4748
    48     MHillas        *fHil;              //! Pointer to MHillas container
    49     MHillasSrc     *fHilSrc;           //! Pointer to MHillasSrc container
    50     MHillasSrc     *fHilAnti;          //! Pointer to MHillasSrc container called MHillasSrcAnti
    51     MHillasExt     *fHilExt;           //! Pointer to MHillasExt container
    52     MParameterD    *fThetaSq;          //! Pointer to MParameterD container called ThetaSq
    53     MParameterD    *fDisp;             //! Pointer to MParameterD container called Disp
    54     MPointingPos   *fPointing;         //! Pointer to MPointingPos container
     49    MHillas        *fHil;               //! Pointer to MHillas container
     50    MHillasSrc     *fHilSrc;            //! Pointer to MHillasSrc container
     51    MHillasSrc     *fHilAnti;           //! Pointer to MHillasSrc container called MHillasSrcAnti
     52    MHillasExt     *fHilExt;            //! Pointer to MHillasExt container
     53    MNewImagePar   *fNewImgPar;         //! Pointer to MHillasExt container
     54    MParameterD    *fThetaSq;           //! Pointer to MParameterD container called ThetaSq
     55    MParameterD    *fDisp;              //! Pointer to MParameterD container called Disp
    5556
    5657    Float_t         fMm2Deg;            //! Conversion factor from mm to deg, from MGeomCam
    … …  
    7576    Double_t GetVal(Int_t i) const;
    7677    TString  GetParam(Int_t i) const;
    77     Double_t GetThetaSq(Double_t c, Double_t wl, Double_t d, Double_t a, MParameterD *par=0) const;
     78    Double_t GetThetaSq(Double_t p, Double_t wl, Double_t d, Double_t a) const;
    7879
    7980public:
Note: See TracChangeset for help on using the changeset viewer.