Ignore:
Timestamp:
07/26/07 12:13:00 (17 years ago)
Author:
tbretz
Message:
*** empty log message ***
File:
1 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}
Note: See TracChangeset for help on using the changeset viewer.