Changeset 7173 for trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc

Timestamp:

06/28/05 10:22:44 (19 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc (modified) (2 diffs)

trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc

@@ -30 +30 @@
 #include "MFMagicCuts.h"
 
+#include <fstream>
+
+#include "MHMatrix.h"
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+#include "MParList.h"
+
+#include "MGeomCam.h"
+#include "MHillasSrc.h"
+#include "MHillasExt.h"
+#include "MParameters.h"
+#include "MPointingPos.h"
+
 ClassImp(MFMagicCuts);
 
 using namespace std;
 
-
 // --------------------------------------------------------------------------
 //
@@ -40 +54 @@
 //
 MFMagicCuts::MFMagicCuts(const char *name, const char *title)
+    : fHil(0), fHilSrc(0), fHilAnti(0), fHilExt(0), fThetaSq(0),
+    fDisp(0), fPointing(0), fMatrix(0), fVariables(30), fThetaCut(kNone),
+    fHadronnessCut(kAll)
 {
     fName  = name  ? name  : "MFMagicCuts";
     fTitle = title ? title : "Class to evaluate the MagicCuts";
-}
+
+    fMatrix = NULL;
+
+    AddToBranchList("MHillas.fWidth");
+    AddToBranchList("MHillas.fLength");
+    AddToBranchList("MHillas.fSize");
+    AddToBranchList("MHillasSrc.fAlpha");
+    AddToBranchList("MHillasSrcAnti.fAlpha");
+    AddToBranchList("MPointingPos.fZd");
+    AddToBranchList("MHillasExt.fMaxDist");
+    AddToBranchList("MHillasExt.fM3Long");
+
+    fVariables[0] =  1.42547;      // Xi
+    fVariables[1] =  0.233773;     // Theta^2
+    fVariables[2] =  0.2539460;    // Area[0]
+    fVariables[3] =  5.2149800;    // Area[1]
+    fVariables[4] =  0.1139130;    // Area[2]
+    fVariables[5] = -0.0889;       // M3long
+    fVariables[6] =  0.18;         // Xi(theta)
+}
+
+// --------------------------------------------------------------------------
+//
+// The theta cut value GetThetaSqCut() is propageted to the parameter list
+// as 'ThetaSquaredCut' as MParameterD so that it can be used somewhere else.
+//
+Int_t MFMagicCuts::PreProcess(MParList *pList)
+{
+    MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
+    if (!cam)
+    {
+        *fLog << err << "MGeomCam not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fMm2Deg = cam->GetConvMm2Deg();
+
+    fThetaSq = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaSquared");
+    if (!fThetaSq)
+        return kFALSE;
+    fDisp = (MParameterD*)pList->FindCreateObj("MParameterD", "Disp");
+    if (!fDisp)
+        return kFALSE;
+
+    // propagate Theta cut to the parameter list
+    // so that it can be used somewhere else
+    MParameterD *thetacut = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaSquaredCut");
+    if (!thetacut)
+        return kFALSE;
+    thetacut->SetVal(GetThetaSqCut());
+    //thetacut->SetReadyToSave();
+
+    MParameterD *m3lcut = (MParameterD*)pList->FindCreateObj("MParameterD", "M3lCut");
+    if (!m3lcut)
+        return kFALSE;
+    m3lcut->SetVal(fVariables[5]);
+
+    MParameterD *dispxi = (MParameterD*)pList->FindCreateObj("MParameterD", "DispXi");
+    if (!dispxi)
+        return kFALSE;
+    dispxi->SetVal(fVariables[0]);
+
+    dispxi = (MParameterD*)pList->FindCreateObj("MParameterD", "DispXiTheta");
+    if (!dispxi)
+        return kFALSE;
+    dispxi->SetVal(fVariables[6]);
+
+    Print();
+
+    if (fMatrix)
+        return kTRUE;
+
+    //-----------------------------------------------------------
+
+    fHil = (MHillas*)pList->FindObject("MHillas");
+    if (!fHil)
+    {
+        *fLog << err << "MHillas not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fHilExt = (MHillasExt*)pList->FindObject("MHillasExt");
+    if (!fHilExt)
+    {
+        *fLog << err << "MHillasExt not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fPointing = (MPointingPos*)pList->FindObject("MPointingPos");
+    if (!fPointing)
+    {
+        *fLog << err << "MPointingPos not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fHilSrc = (MHillasSrc*)pList->FindObject("MHillasSrc");
+    if (!fHilSrc)
+    {
+        *fLog << err << "MHillasSrc not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    if (fThetaCut&kOff)
+    {
+        fHilAnti = (MHillasSrc*)pList->FindObject("MHillasSrcAnti", "MHillasSrc");
+        if (!fHilAnti)
+        {
+            *fLog << warn << "MHillasSrcAnti [MHillasSrc] not found... aborting." << endl;
+            return kFALSE;
+        }
+    }
+
+    return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+// Returns the mapped value from the Matrix
+//
+Double_t MFMagicCuts::GetVal(Int_t i) const
+{
+    return (*fMatrix)[fMap[i]];
+}
+
+// --------------------------------------------------------------------------
+//
+// You can use this function if you want to use a MHMatrix instead of the
+// given containers. This function adds all necessary columns to the
+// given matrix. Afterward you should fill the matrix with the corresponding
+// data (eg from a file by using MHMatrix::Fill). If you now loop
+// through the matrix (eg using MMatrixLoop) MEnergyEstParam::Process
+// will take the values from the matrix instead of the containers.
+//
+void MFMagicCuts::InitMapping(MHMatrix *mat)
+{
+    if (fMatrix)
+      return;
+
+    fMatrix = mat;
+
+    fMap[kESize]   = fMatrix->AddColumn("log10(MHillas.fSize)");
+    fMap[kEArea]   = fMatrix->AddColumn("MHillas.GetArea*MGeomCam.fConvMm2Deg*MGeomCam.fConvMm2Deg");
+    fMap[kEM3Long] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");
+    fMap[kEWdivL]  = fMatrix->AddColumn("MHillas.fWidth/MHillas.fLength");
+    fMap[kEZd]     = fMatrix->AddColumn("MPointingPos.GetZdRad");
+
+    fMap[kEAlpha]  = fMatrix->AddColumn("MHillasSrc.fAlpha");
+    fMap[kEDist]   = fMatrix->AddColumn("MHillasSrc.fDist*MGeomCam.fConvMm2Deg");
+    //fMap[kDCA]    = fMatrix->AddColumn("MHillasSrc.fDCA*MGeomCam.fConvMm2Deg");
+    if (fThetaCut&kOff)
+    {
+        fMap[kEAlphaAnti]  = fMatrix->AddColumn("MHillasSrcAnti.fAlpha");
+        fMap[kEDistAnti]   = fMatrix->AddColumn("MHillasSrcAnti.fDist*MGeomCam.fConvMm2Deg");
+        fMap[kEM3LongAnti] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrcAnti.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");
+        //fMap[kDCAAnti]   = fMatrix->AddColumn("MHillasSrcAnti.fDCA*MGeomCam.fConvMm2Deg");
+    }
+
+    //fMap[kLength] = fMatrix->AddColumn("MHillas.fLength*MGeomCam.fConvMm2Deg");
+    //fMap[kWidth]  = fMatrix->AddColumn("MHillas.fWidth*MGeomCam.fConvMm2Deg");
+}
+
+// --------------------------------------------------------------------------
+//
+// Returns theta squared based on the formula:
+//    p := c*(1-width/length)
+//    return d*d + p*p - 2*d*p*cos(a*TMath::DegToRad());
+//
+// If par!=NULL p is stored in this MParameterD
+//
+Double_t MFMagicCuts::GetThetaSq(Double_t c, Double_t wl, Double_t d, Double_t a, MParameterD *par) const
+{
+    // wl = width/l [1]
+    // d  = dist    [deg]
+    // a  = alpha   [deg]
+    const Double_t p = c*(1-wl);
+    if (par)
+        par->SetVal(p);
+    return d*d + p*p - 2*d*p*cos(a*TMath::DegToRad());
+}
+
+// --------------------------------------------------------------------------
+//
+// Returns squared cut value in theta: fVariables[1]^2
+//
+Double_t MFMagicCuts::GetThetaSqCut() const
+{
+    return fVariables[1]*fVariables[1];
+}
+
+// ---------------------------------------------------------------------------
+//
+// Evaluate dynamical magic-cuts
+//
+Int_t MFMagicCuts::Process()
+{
+    // Get some variables
+    const Double_t wdivl  = fMatrix ? GetVal(kEWdivL) : fHil->GetWidth()/fHil->GetLength();
+    const Double_t lgsize = fMatrix ? GetVal(kESize)  : TMath::Log10(fHil->GetSize());
+    const Double_t zd     = fMatrix ? GetVal(kEZd)    : fPointing->GetZdRad();
+
+    // For simplicity
+    const Double_t *c = fVariables.GetArray();
+
+    // Value for Theta cut
+    const Double_t cut = GetThetaSqCut(); //c[1]*c[1];
+    const Double_t xi  = (c[0]+c[8]*(lgsize-c[7])*(lgsize-c[7])) - c[6]*zd*zd;
+
+    // Default if we return before the end
+    fResult = kFALSE;
+
+    // ------------------- Most efficient cut -----------------------
+    // ---------------------- Theta cut ON --------------------------
+    const Double_t dist    = fMatrix ? GetVal(kEDist)  : fHilSrc->GetDist()*fMm2Deg;
+    const Double_t alpha   = fMatrix ? GetVal(kEAlpha) : fHilSrc->GetAlpha();
+
+    const Double_t thetasq = GetThetaSq(xi, wdivl, dist, alpha, fDisp);
+
+    fThetaSq->SetVal(thetasq);
+
+    if (fThetaCut&kOn && thetasq >= cut)
+        return kTRUE;
+
+    // --------------------- Efficient  cut -------------------------
+    // ---------------------- Area/M3l cut --------------------------
+    if (fHadronnessCut&kArea)
+    {
+        const Double_t area = fMatrix ? GetVal(kEArea) : fHil->GetArea()*fMm2Deg*fMm2Deg;
+        const Double_t A    = lgsize>c[3] ? c[2] : c[2]*(1 - c[4]*(lgsize-c[3])*(lgsize-c[3]));
+        if (area>=A)
+            return kTRUE;
+    }
+    if (fHadronnessCut&kM3Long)
+    {
+        const Double_t m3long = fMatrix ? GetVal(kEM3Long) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha());
+        if (m3long<=c[5])
+            return kTRUE;
+    }
+
+    // ------------------- Least efficient cut -----------------------
+    // ---------------------- Theta cut OFF --------------------------
+    if (fThetaCut&kOff)
+    {
+        const Double_t distanti    = fMatrix ? GetVal(kEDistAnti)   : fHilAnti->GetDist()*fMm2Deg;
+        const Double_t alphaanti   = fMatrix ? GetVal(kEAlphaAnti)  : fHilAnti->GetAlpha();
+        const Double_t thetasqanti = GetThetaSq(xi, wdivl, distanti, alphaanti);
+        const Double_t m3lanti     = fMatrix ? GetVal(kEM3LongAnti) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha());
+
+        if (thetasqanti<cut && (fHadronnessCut&kM3Long || m3lanti>c[5]))
+            return kTRUE;
+    }
+
+    fResult = kTRUE;
+
+    return kTRUE;
+}
+
+void MFMagicCuts::SetVariables(const TArrayD &arr)
+{
+    fVariables = arr;
+}
+
+TString MFMagicCuts::GetParam(Int_t i) const
+{
+    if (i>=fVariables.GetSize())
+        return "";
+
+    return Form("Param[%2d]=%+f", i, fVariables[i]);
+}
+
+void MFMagicCuts::Print(Option_t *o) const
+{
+    *fLog << all << GetDescriptor() << ":" << setiosflags(ios::left) << endl;
+
+    *fLog << "Using Theta cut: ";
+    switch (fThetaCut)
+    {
+    case kNone:
+        *fLog << "none" << endl;
+        break;
+    case kOn:
+        *fLog << "on" << endl;
+        break;
+    case kOff:
+        *fLog << "off" << endl;
+        break;
+    case kWobble:
+        *fLog << "on and off (wobble)" << endl;
+        break;
+    }
+    *fLog << "Using hadronness cut: ";
+    switch (fHadronnessCut)
+    {
+    case kNoCut:
+        *fLog << "none" << endl;
+        break;
+    case kArea:
+        *fLog << "area" << endl;
+        break;
+    case kM3Long:
+        *fLog << "m3long" << endl;
+        break;
+    case kAll:
+        *fLog << "all" << endl;
+        break;
+    }
+    *fLog << "Filter is" << (IsInverted()?"":" not") << " inverted." << endl;
+
+    const Int_t n = fVariables.GetSize();
+    for (int i=0; i<n; i+=3)
+    {
+        *fLog << setw(25) << GetParam(i);
+        *fLog << setw(25) << GetParam(i+1);
+        *fLog << setw(25) << GetParam(i+2);
+        *fLog << endl;
+    }
+    *fLog << setiosflags(ios::right);
+}
+
+Bool_t MFMagicCuts::CoefficentsRead(const char *fname)
+{
+    ifstream fin(fname);
+    if (!fin)
+    {
+        *fLog << err << "Cannot open file " << fname << ": ";
+        *fLog << strerror(errno) << endl;
+        return kFALSE;
+    }
+
+    for (int i=0; i<fVariables.GetSize(); i++)
+    {
+        fin >> fVariables[i];
+        if (!fin)
+        {
+            *fLog << err << "ERROR - Not enough coefficients in file " << fname << endl;
+            return kERROR;
+        }
+    }
+    return kTRUE;
+}
+
+Bool_t MFMagicCuts::CoefficentsWrite(const char *fname) const
+{
+    ofstream fout(fname);
+    if (!fout)
+    {
+        *fLog << err << "Cannot open file " << fname << ": ";
+        *fLog << strerror(errno) << endl;
+        return kFALSE;
+    }
+
+    for (int i=0; i<fVariables.GetSize(); i++)
+        fout << setw(11) << fVariables[i] << endl;
+
+    return kTRUE;
+}
+
+Int_t MFMagicCuts::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
+{
+    if (MFilter::ReadEnv(env, prefix, print)==kERROR)
+        return kERROR;
+
+
+    Bool_t rc = kFALSE;
+    if (IsEnvDefined(env, prefix, "ThetaCut", print))
+    {
+        TString str = GetEnvValue(env, prefix, "ThetaCut", "");
+        str.ToLower();
+        str=str.Strip(TString::kBoth);
+
+        if (str==(TString)"nocut" || str==(TString)"none")
+            fThetaCut = kNone;
+        if (str==(TString)"on")
+            fThetaCut = kOn;
+        if (str==(TString)"off")
+            fThetaCut = kOff;
+        if (str==(TString)"wobble")
+            fThetaCut = kWobble;
+        rc = kTRUE;
+    }
+    if (IsEnvDefined(env, prefix, "HadronnessCut", print))
+    {
+        TString str = GetEnvValue(env, prefix, "HadronnessCut", "");
+        str.ToLower();
+        str=str.Strip(TString::kBoth);
+
+        if (str==(TString)"nocut" || str==(TString)"none")
+            fHadronnessCut = kNoCut;
+        if (str==(TString)"area")
+            fHadronnessCut = kArea;
+        if (str==(TString)"m3long")
+            fHadronnessCut = kM3Long;
+        if (str==(TString)"all")
+            fHadronnessCut = kAll;
+
+        rc = kTRUE;
+    }
+
+    if (IsEnvDefined(env, prefix, "File", print))
+    {
+        const TString fname = GetEnvValue(env, prefix, "File", "");
+        if (!CoefficentsRead(fname))
+            return kERROR;
+
+        return kTRUE;
+    }
+
+    for (int i=0; i<fVariables.GetSize(); i++)
+    {
+        if (IsEnvDefined(env, prefix, Form("Param%d", i), print))
+        {
+            // It is important that the last argument is a floating point
+            fVariables[i] = GetEnvValue(env, prefix, Form("Param%d", i), 0.0);
+            rc = kTRUE;
+        }
+    }
+    return rc;
+    //return kTRUE; // means: can use default values
+    //return rc;  // means: require something in resource file
+}