Ignore:
Timestamp:
07/26/07 12:13:00 (17 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

    r8633 r8646  
    5151//
    5252//    c[0],  c[5], c[6], c[7],  c[8], c[9]:
    53 //          xi          = c[0]+c[6]*pow(leakage1, c[7]);
     53//          xi          = c[0] + c[8]*slope + c[9]*leak +
     54//                        (lgsize>c[10])*c[11]*(lgsize-c[10])^2;
    5455//          p           = xi*(1-width/length);
    5556//          sign1       = m3long-c[5]
    56 //          sign2       = (dist-c[9])*c[8]-slope
     57//          sign2       = (dist-c[7])*c[6]-slope
    5758//          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
    6159//          thetasq     = disp^2 + dist^2 - 2*disp*dist*alpha
    62 //          antithetasq = antidisp^2 + antidist^2 - 2*antidisp*antidist*antialpha
     60//
     61//    And the values with respect to the antisource position respectively.
     62//
    6363//
    6464//    c[1]:
    … …  
    7171//
    7272//    HadronnessCut:
    73 //    c[10], c[11]:
    74 //          had       <= c[10]
    75 //          10^lgsize >= c[11]
     73//    c[13], c[14]:
     74//          had       <= c[13]
     75//          10^lgsize >= c[14]
    7676//
    7777//
    … …  
    8181//   HadronnessCut:
    8282//    - none/nocut: No area cut
    83 //    - area:       Area cut
    84 //    - all:        same as area
     83//    - area:       Area cut <default>
     84//    - hadronness: Cut in size and hadronness (c[10], c[11])
     85//    - all:        area + hadronness
    8586//
    8687//   ThetaCut:
    87 //    - none/nocut: no theta cut
     88//    - none/nocut: no theta cut <default>
    8889//    - on:         cut in theta only
    8990//    - off:        anti-cut in anti-theta only
    9091//    - wobble:     same as on|off (both)
    9192//
     93//   CalcDisp:
     94//    - yes:        Disp is calculated as defined above <default>
     95//    - no:         abs(Disp.fVal) from the parameter list is used instead
     96//
     97//   CalcGhostbuster:
     98//    - yes:        The ghostbuster is calculated as defined above <default>
     99//    - no:         Ghostbuster.fVal<c[12] is used as ghostbusting condition
    92100//
    93101//  Input Container:
    … …  
    99107//   MHillasSrc
    100108//   [MHillasSrcAnti [MHillasSrc]]
     109//   [Disp [MParameterD]]
     110//   [Ghostbuster [MParameterD]]
    101111//
    102112//  Output:
    … …  
    138148MFMagicCuts::MFMagicCuts(const char *name, const char *title)
    139149    : fHil(0), fHilSrc(0), fHilAnti(0), fHilExt(0), fNewImgPar(0),
    140     fThetaSq(0), fDisp(0), fHadronness(0), fMatrix(0), fVariables(30),
    141     fThetaCut(kNone), fHadronnessCut(kArea)
     150    fThetaSq(0), fDisp(0), fHadronness(0), fMatrix(0), fVariables(20),
     151    fThetaCut(kNone), fHadronnessCut(kArea), fCalcDisp(kTRUE),
     152    fCalcGhostbuster(kTRUE)
    142153{
    143154    fName  = name  ? name  : "MFMagicCuts";
    … …  
    159170    AddToBranchList("MNewImagePar.fLeakage1");
    160171    AddToBranchList("Hadronness.fVal");
    161 /*
    162     fVariables[0] =  1.42547;      // Xi
    163     fVariables[1] =  0.233773;     // Theta^2
    164     fVariables[2] =  0.2539460;    // Area[0]
    165     fVariables[3] =  5.2149800;    // Area[1]
    166     fVariables[4] =  0.1139130;    // Area[2]
    167     fVariables[5] = -0.0889;       // M3long
    168     fVariables[6] =  0.18;         // Xi(theta)
    169     fVariables[7] =  0.18;         // Xi(theta)
    170     */
     172    AddToBranchList("Disp.fVal");
     173    AddToBranchList("Ghostbuster.fVal");
    171174}
    172175
    … …  
    190193    if (!fThetaSq)
    191194        return kFALSE;
    192     fDisp = (MParameterD*)pList->FindCreateObj("MParameterD", "Disp");
     195
     196    if (!fCalcDisp)
     197        fDisp = (MParameterD*)pList->FindObject("Disp", "MParameterD");
     198    else
     199        fDisp = (MParameterD*)pList->FindCreateObj("MParameterD", "Disp");
    193200    if (!fDisp)
    194         return kFALSE;
     201    {
     202        *fLog << err << "Disp [MParameterD] not found... aborting." << endl;
     203        return kFALSE;
     204    }
     205
     206    if (!fCalcGhostbuster)
     207        fGhostbuster = (MParameterD*)pList->FindObject("Ghostbuster", "MParameterD");
     208    else
     209        fGhostbuster = (MParameterD*)pList->FindCreateObj("MParameterD", "Ghostbuster");
     210    if (!fGhostbuster)
     211    {
     212        *fLog << err << "Ghostbuster [MParameterD] not found... aborting." << endl;
     213        return kFALSE;
     214    }
    195215
    196216    // propagate Theta cut to the parameter list
    … …  
    297317    fMap[kEM3Long]  = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");
    298318
    299     fMap[kESrcSign] = fMatrix->AddColumn("sign(MHillasSrc.fCosDeltaAlpha)");
     319    fMap[kESign]    = fMatrix->AddColumn("sign(MHillasSrc.fCosDeltaAlpha)");
    300320
    301321    fMap[kESlope]   = fMatrix->AddColumn("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/MGeomCam.fConvMm2Deg");
     322
     323    if (!fCalcDisp)
     324        fMap[kEDisp] = fMatrix->AddColumn("abs(Disp.fVal)");
     325
     326    if (!fCalcGhostbuster)
     327        fMap[kEGhostbuster] = fMatrix->AddColumn("Ghostbuster.fVal");
    302328
    303329    if (fThetaCut&kOff)
    … …  
    321347// If par!=NULL p is stored in this MParameterD
    322348//
    323 Double_t MFMagicCuts::GetThetaSq(Double_t p, Double_t wl, Double_t d, Double_t a) const
     349Double_t MFMagicCuts::GetThetaSq(Double_t p, Double_t d, Double_t a) const
    324350{
    325351    // wl = width/l [1]
    … …  
    338364}
    339365
    340 // ---------------------------------------------------------------------------
    341 //
    342 // Evaluate dynamical magic-cuts
    343 //
    344 Int_t MFMagicCuts::Process()
    345 {
     366// --------------------------------------------------------------------------
     367//
     368// Return abs(Disp.fVal) if disp calculation is switched off.
     369// Otherwise return (c0+c6*leak^c7)*(1-width/length)
     370//
     371Double_t MFMagicCuts::GetDisp(Double_t slope, Double_t lgsize) const
     372{
     373    if (!fCalcDisp)
     374        return fMatrix ? GetVal(kEDisp) : TMath::Abs(fDisp->GetVal());
     375
    346376    // Get some variables
    347     const Double_t wdivl  = fMatrix ? GetVal(kEWdivL)   : fHil->GetWidth()/fHil->GetLength();
    348     const Double_t lgsize = fMatrix ? GetVal(kESize)    : TMath::Log10(fHil->GetSize());
    349     const Double_t leak   = fMatrix ? GetVal(kELeakage) : TMath::Log10(fNewImgPar->GetLeakage1()+1);
     377    const Double_t wdivl = fMatrix ? GetVal(kEWdivL)   : fHil->GetWidth()/fHil->GetLength();
     378    const Double_t leak  = fMatrix ? GetVal(kELeakage) : TMath::Log10(fNewImgPar->GetLeakage1()+1);
    350379
    351380    // For simplicity
    352381    const Double_t *c = fVariables.GetArray();
    353382
     383    // As rule for root or MDataPhrase:
     384    //   ((M[3]>[3])*[4]*(M[3]-[3])^2 + [2]*M[2] + [1]*M[1] + [0])*M[0]
     385    //
     386    Double_t xi = c[0] + c[8]*slope + c[9]*leak;
     387    if (lgsize>c[10])
     388        xi += (lgsize-c[10])*(lgsize-c[10])*c[11];
     389
     390    const Double_t disp = xi*(1-wdivl);
     391
     392    return disp;
     393}
     394
     395Bool_t MFMagicCuts::IsGhost(Double_t m3long, Double_t slope, Double_t dist) const
     396{
     397    // For simplicity
     398    const Double_t *c = fVariables.GetArray();
     399
     400    if (!fCalcGhostbuster)
     401        return (fMatrix ? GetVal(kEGhostbuster) : fGhostbuster->GetVal())<c[12];
     402
     403    // Do Ghostbusting with slope and m3l
     404    const Bool_t sign1 = m3long < c[5];
     405    const Bool_t sign2 = slope  > (dist-c[7])*c[6];
     406
     407    return sign1 || sign2;
     408}
     409
     410// ---------------------------------------------------------------------------
     411//
     412// Evaluate dynamical magic-cuts
     413//
     414Int_t MFMagicCuts::Process()
     415{
     416    // For simplicity
     417    const Double_t *c = fVariables.GetArray();
     418
     419    // Default if we return before the end
     420    fResult = kFALSE;
     421
    354422    // Value for Theta cut (Disp parametrization)
    355423    const Double_t cut  = GetThetaSqCut();
    356     const Double_t xi   = c[0]+c[6]*pow(leak, c[7]);
    357     const Double_t disp = xi*(1-wdivl);
    358 
    359     // Default if we return before the end
    360     fResult = kFALSE;
    361424
    362425    // ------------------- Most efficient cut -----------------------
    363426    // ---------------------- Theta cut ON --------------------------
    364     const Double_t dist   = fMatrix ? GetVal(kEDist)    : fHilSrc->GetDist()*fMm2Deg;
    365     const Double_t alpha  = fMatrix ? GetVal(kEAlpha)   : fHilSrc->GetAlpha();
    366     const Double_t sign   = fMatrix ? GetVal(kESrcSign) : MMath::Sgn(fHilSrc->GetCosDeltaAlpha());
    367     const Double_t m3long = fMatrix ? GetVal(kEM3Long)  : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, 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];
     427    const Double_t dist   = fMatrix ? GetVal(kEDist)   : fHilSrc->GetDist()*fMm2Deg;
     428    const Double_t alpha  = fMatrix ? GetVal(kEAlpha)  : fHilSrc->GetAlpha();
     429    const Double_t m3long = fMatrix ? GetVal(kEM3Long) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha());
     430    const Double_t slope  = fMatrix ? GetVal(kESlope)  : fHilExt->GetSlopeLong()/TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha());
     431    const Double_t sign   = fMatrix ? GetVal(kESign)   : MMath::Sgn(fHilSrc->GetCosDeltaAlpha());
     432    const Double_t lgsize = fMatrix ? GetVal(kESize)   : TMath::Log10(fHil->GetSize());
    373433
    374434    // Calculate disp including sign
    375     const Double_t p = sign1<0 || sign2<0 ? -disp : disp;
    376 
    377     // Align disp along source direction depending on third moment
    378     //const Double_t p = TMath::Sign(disp, m3long-c[5]);
     435    const Double_t disp  = GetDisp(slope, lgsize);
     436    const Double_t ghost = IsGhost(m3long, slope, dist);
     437
     438    const Double_t p = ghost ? -disp : disp;
    379439
    380440    // Align sign of disp along shower axis like m3long
    … …  
    382442
    383443    // Calculate corresponding Theta^2
    384     const Double_t thetasq = GetThetaSq(p, wdivl, dist, alpha);
     444    const Double_t thetasq = GetThetaSq(p, dist, alpha);
    385445    fThetaSq->SetVal(thetasq);
    386446
    … …  
    405465    {
    406466        const Double_t had = fMatrix ? GetVal(kEHadronness) : fHadronness->GetVal();
    407         if (had>c[10])
     467        if (had>c[13])
    408468            return kTRUE;
    409469
    410         if (TMath::Power(10, lgsize)<c[11])
     470        if (TMath::Power(10, lgsize)<c[14])
    411471            return kTRUE;
    412472    }
    … …  
    422482        const Double_t slopeanti   = fMatrix ? GetVal(kESlopeAnti)  : fHilExt->GetSlopeLong()/TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha());
    423483
    424         // Do Ghostbusting with slope and m3l
    425         const Double_t sign3 = (distanti-c[9])*c[8]-slopeanti;
    426         const Double_t sign4 = m3lanti-c[5];
    427         const Double_t panti = sign3<0 || sign4<0 ? -disp : disp;
     484        const Double_t dispanti  = GetDisp(slopeanti, lgsize);
     485        const Double_t ghostanti = IsGhost(m3lanti, slopeanti, lgsize);
     486
     487        const Double_t panti = ghostanti ? -dispanti : dispanti;
    428488
    429489        // Align disp along source direction depending on third moment
    430         //const Double_t panti       = TMath::Sign(disp, m3lanti-c[5]);
    431         const Double_t thetasqanti = GetThetaSq(panti, wdivl, distanti, alphaanti);
     490        const Double_t thetasqanti = GetThetaSq(panti, distanti, alphaanti);
    432491
    433492        if (thetasqanti<cut)
    … …  
    494553        break;
    495554    }
     555    if (fCalcDisp)
     556        *fLog << "Disp is calculated from c0, c7, c8." << endl;
     557    else
     558        *fLog << "Disp.fVal from the parameter list is used as disp." << endl;
     559    if (fCalcGhostbuster)
     560        *fLog << "Ghostbusting is calculated from c5, c6, c7." << endl;
     561    else
     562        *fLog << "Ghostbuster.fVal from the parameter list is used for ghostbusting." << endl;
     563
    496564    *fLog << "Filter is" << (IsInverted()?"":" not") << " inverted." << endl;
    497565
    … …  
    568636        rc = kTRUE;
    569637    }
     638
    570639    if (IsEnvDefined(env, prefix, "HadronnessCut", print))
    571640    {
    … …  
    586655    }
    587656
     657    if (IsEnvDefined(env, prefix, "CalcDisp", print))
     658    {
     659        fCalcDisp = GetEnvValue(env, prefix, "CalcDisp", fCalcDisp);
     660        rc = kTRUE;
     661    }
     662
     663    if (IsEnvDefined(env, prefix, "CalcGhostbuster", print))
     664    {
     665        fCalcGhostbuster = GetEnvValue(env, prefix, "CalcGhostbuster", fCalcGhostbuster);
     666        rc = kTRUE;
     667    }
     668
    588669    if (IsEnvDefined(env, prefix, "File", print))
    589670    {
    … …  
    605686    }
    606687    return rc;
    607     //return kTRUE; // means: can use default values
    608     //return rc;  // means: require something in resource file
    609 }
     688}

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

    r8619 r8646  
    4343    enum {
    4444        kESize, kEAlpha, kEAlphaAnti, kEArea, kEDist, kEDistAnti,
    45         kEM3Long, kEM3LongAnti, kEWdivL, kELeakage, kESrcSign,
    46         kESlope, kESlopeAnti, kEHadronness,
     45        kEM3Long, kEM3LongAnti, kEWdivL, kELeakage, kESlope,
     46        kESlopeAnti, kEHadronness, kESign, kEDisp, kEGhostbuster,
    4747        kLastElement
    4848    };
    … …  
    5555    MParameterD    *fThetaSq;           //! Pointer to MParameterD container called ThetaSq
    5656    MParameterD    *fDisp;              //! Pointer to MParameterD container called Disp
     57    MParameterD    *fGhostbuster;       //! Pointer to MParameterD container called Ghostbuster
    5758    MParameterD    *fHadronness;        //! Pointer to MParameterD container called Hadronness
    5859
    … …  
    6768    ThetaCut_t      fThetaCut;          // Which kind of theta cut should be evaluated
    6869    HadronnessCut_t fHadronnessCut;     // Which kind of hadronness cut should be evaluated
     70    Bool_t          fCalcDisp;          // Should we use Disp from the parameterlist?
     71    Bool_t          fCalcGhostbuster;   // Should we use Ghostbuster from the parameterlist?
    6972
    7073    // MTask
    … …  
    7679
    7780    // MFMagicCuts
     81    Double_t GetDisp(Double_t slope, Double_t lgsize) const;
     82    Bool_t   IsGhost(Double_t m3long, Double_t slope, Double_t dist) const;
    7883    Double_t GetVal(Int_t i) const;
    7984    TString  GetParam(Int_t i) const;
    80     Double_t GetThetaSq(Double_t p, Double_t wl, Double_t d, Double_t a) const;
     85    Double_t GetThetaSq(Double_t p, Double_t d, Double_t a) const;
    8186
    8287public:
    … …  
    8994    void   SetThetaCut(ThetaCut_t c) { fThetaCut=c; }
    9095    void   SetHadronnessCut(HadronnessCut_t c) { fHadronnessCut=c; }
     96    void   SetCalcDisp(Bool_t b=kTRUE) { fCalcDisp=b; }
     97    void   SetCalcGhostbuster(Bool_t b=kTRUE) { fCalcGhostbuster=b; }
    9198
    9299    // MFMagicCuts
Note: See TracChangeset for help on using the changeset viewer.