source: trunk/MagicSoft/Mars/mfilter/MFSoftwareTrigger.cc@ 6802

/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Wolfgang Wittek, 04/2003 <mailto:wittek@mppmu.mpg.de>
!   Author(s): Thomas Bretz, 04/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2005
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MFSoftwareTrigger
//
//  This is a class to evaluate a software trigger
//
//  The number of required pixels is setup by:
//   SetNumNeighbors(4)   // default
//
//  The Threshold is setup by:
//    SetThreshold(5)     // default
//
//  Time window for coincidence:
//    SetTimeWindow(3.3)
// To switch off time-coincidence use
//    SetTimeWindow(-1)
//   or
//    SetTimeWindow()
//
//  kSinglePixelsNeighbors <default>
//  ----------------------
//    Checks whether there is at least one pixel above threshold which
//    has fNumNeighbors direct neighbors which are also above threshold.
//    The outermost ring of pixels is ignord. Only 'used' pixels are
//    taken into account.
//
//  kAnyPattern 
//  -----------
//    Checks whether it find a cluster of pixels above fThreshold which
//    has a size bigger/equal than fNumNeighbors. The layout (pattern) of
//    the cluster is ignored. Unmapped pixels are ignored.
//
//  WARNING: Using trigger type kAllPattern resets the BIT(21) bit
//           of all pixels in MCerPhotEvt
//
//
//  Input:
//    MCerPhotEvt
//    MArrivalTime
//    MGeomCam
//
/////////////////////////////////////////////////////////////////////////////
#include "MFSoftwareTrigger.h"

#include "MLog.h"
#include "MLogManip.h"

#include "MParList.h"

#include "MGeomPix.h"
#include "MGeomCam.h"

#include "MCerPhotEvt.h"
#include "MArrivalTime.h"

ClassImp(MFSoftwareTrigger);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor.
//
MFSoftwareTrigger::MFSoftwareTrigger(const char *name, const char *title)
    : fCam(NULL), fEvt(NULL), fTme(NULL), fThreshold(5),
    fTimeWindow(0.5), fNumNeighbors(4), fType(kSinglePixelNeighbors)
{
    fName  = name  ? name  : "MFSoftwareTrigger";
    fTitle = title ? title : "Filter for software trigger";
}

// --------------------------------------------------------------------------
//
// This function recursively finds all pixels of one island and sets
// BIT(14) for the pixel.
//
//  1) Check whether a pixel with the index idx exists, is unused
//     and has not yet BIT(14) set
//  2) Set BIT(14) to the pixel
//  3) Loop over all its neighbors taken from the geometry geom. For all
//     neighbors recursively call this function (CountPixels)
//  4) Sum the number of valid neighbors newly found
//
// Returns the size of the cluster
//
Int_t MFSoftwareTrigger::CountPixels(Int_t idx, Float_t tm0) const
{
    // Try to get the pixel information of a pixel with this index
    MCerPhotPix *pix = fEvt->GetPixById(idx);

    // If a pixel with this index is not existing... do nothing.
    if (!pix)
        return 0;

    // If pixel already assigned to a cluster
    if (pix->TestBit(kWasChecked))
        return 0;

    if (pix->IsPixelUnmapped())
        return 0;

    // Assign the new island number num to this used pixel
    pix->SetBit(kWasChecked);

    // Get the size of this pixel and check threshold
    const Double_t size = pix->GetNumPhotons()*fCam->GetPixRatio(idx);
    if (size<fThreshold)
        return 0;

    const Float_t tm1 = fTme ? (*fTme)[idx] : 0;
    if (fTme && TMath::Abs(tm1-tm0)>fTimeWindow)
        return 0;

    //pix->SetBit(kAboveThreshold);

    Int_t num = 1;

    // Get the geometry information (neighbors) of this pixel
    const MGeomPix &gpix = (*fCam)[idx];

    // Now do the same with all its neighbors and sum the
    // sizes which they correspond to
    const Int_t n = gpix.GetNumNeighbors();
    for (int i=0; i<n; i++)
        num += CountPixels(gpix.GetNeighbor(i), tm1);

    // return size of this (sub)cluster
    return num;
}
/*
Int_t MFSoftwareTrigger::CountCoincidencePixels(Int_t idx) const
{
    // Try to get the pixel information of a pixel with this index
    MCerPhotPix *pix = fEvt->GetPixById(idx);

    // If a pixel with this index is not existing... do nothing.
    if (!pix)
        return 0;

    // If pixel already assigned to a cluster
    if (pix->TestBit(kWasChecked))
        return 0;

    if (pix->IsPixelUnmapped())
        return 0;

    // Assign the new island number num to this used pixel
    pix->SetBit(kWasChecked);

    // Get the size of this pixel and check threshold
    const Double_t size = pix->GetNumPhotons();
    if (size<fThreshold)
        return 0;

    Int_t num = 1;

    // Get the geometry information (neighbors) of this pixel
    const MGeomPix &gpix = (*fCam)[idx];

    // Now do the same with all its neighbors and sum the
    // sizes which they correspond to
    const Int_t n = gpix.GetNumNeighbors();
    for (int i=0; i<n; i++)
        num += CountPixels(gpix.GetNeighbor(i));

    // return size of this (sub)cluster
    return num;
}
*/
void MFSoftwareTrigger::ResetBits(Int_t bits) const
{
    TObject *obj=0;

    TIter Next(*fEvt);
    while ((obj=Next()))
        obj->ResetBit(bits);
}

// --------------------------------------------------------------------------
//
// Check if there is at least one pixel which fullfills the condition
//
Bool_t MFSoftwareTrigger::ClusterTrigger() const
{
    ResetBits(kWasChecked);

    // Reset bit
    MCerPhotPix *pix=0;

    // We could loop over all indices which looks more straight
    // forward but should be a lot slower (assuming zero supression)
    TIter Next(*fEvt);
    while ((pix=static_cast<MCerPhotPix*>(Next())))
    {
        // Check if trigger condition is fullfilled for this pixel
        const Int_t idx = pix->GetPixId();
        if (CountPixels(idx, (*fTme)[idx]) >= fNumNeighbors)
            return kTRUE;
    }

    return kFALSE;
}
/*
Int_t MFSoftwareTrigger::CheckCoincidence(Int_t idx, Float_t tm0) const
{
    // Try to get the pixel information of a pixel with this index
    MCerPhotPix *pix = fEvt->GetPixById(idx);

    // If a pixel with this index is not existing... do nothing.
    if (!pix)
        return 0;

    // If pixel already assigned to a cluster
    if (pix->TestBit(kWasChecked))
        return 0;

    if (pix->IsPixelUnmapped())
        return 0;

    // Assign the new island number num to this used pixel
    pix->SetBit(kWasChecked);

    const Double_t size = pix->GetNumPhotons();
    if (size<fThreshold)
        return 0;

    const Float_t tm1 = (*fTme)[idx];
    if (TMath::Abs(tm1-tm0)>fTimeWindow)
        return 0;

    pix->SetBit(kIsCoincident);


    Int_t num = 1;

    // Get the geometry information (neighbors) of this pixel
    const MGeomPix &gpix = (*fCam)[idx];

    Int_t cnt = 0;

    // Now do the same with all its neighbors and sum the
    // sizes which they correspond to
    const Int_t n = gpix.GetNumNeighbors();
    for (int i=0; i<n; i++)
    {
        const Int_t rc = CheckCoincidence(gpix.GetNeighbor(i), tm0);
        if (fEvt->GetPixById(gpix.GetNeighbor(i))->TestBit(kIsCoincident))
            cnt++;
        num += rc;
    }

    // return size of this (sub)cluster
    return cnt<2 ? 0 : num;

}

Bool_t MFSoftwareTrigger::MagicLvl1Trigger() const
{
    // Reset bit
    MCerPhotPix *pix=0;

    // We could loop over all indices which looks more straight
    // forward but should be a lot slower (assuming zero supression)
    TIter Next(*fEvt);
    while ((pix=static_cast<MCerPhotPix*>(Next())))
    {
        ResetBits(kWasChecked|kIsCoincident);

        const Int_t idx = pix->GetPixId();
        if (CheckCoincidence(idx, (*fTme)[idx])<fNumNeighbors)
            continue;

        return kTRUE;
    }
    return kFALSE;
}
*/

Bool_t MFSoftwareTrigger::CheckPixel(const MCerPhotPix &pix) const
{
    if (!pix.IsPixelUsed())
        return kFALSE;

    if (pix.GetNumPhotons()*fCam->GetPixRatio(pix.GetPixId())<fThreshold)
        return kFALSE;

    if ((*fCam)[pix.GetPixId()].IsInOutermostRing())
        return kFALSE;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Software single pixel coincidence trigger
//
Bool_t MFSoftwareTrigger::SwTrigger() const
{
    const Int_t entries = fEvt->GetNumPixels();
 
    for (Int_t i=0; i<entries; i++)
    {
        const MCerPhotPix &pix0 = (*fEvt)[i];
        if (!CheckPixel(pix0))
            continue;

        Int_t num = 1;

        const MGeomPix &gpix = (*fCam)[pix0.GetPixId()];

        const Int_t nneighbors = gpix.GetNumNeighbors();
        for (Int_t n=0; n<nneighbors; n++)
        {
            const Int_t idx1 = gpix.GetNeighbor(n);
            if (!CheckPixel(*fEvt->GetPixById(idx1)))
                continue;

            if (fTme)
            {
                const Float_t t0 = (*fTme)[pix0.GetPixId()];
                const Float_t t1 = (*fTme)[idx1];

                if (TMath::Abs(t0-t1)>fTimeWindow)
                    continue;
            }

            if (++num==fNumNeighbors)
                return kTRUE;
        }
    }
    return kFALSE;
}

// --------------------------------------------------------------------------
//
// Request pointer to MCerPhotEvt and MGeomCam from paremeter list
//
Int_t MFSoftwareTrigger::PreProcess(MParList *pList)
{
    fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
    if (!fEvt)
    {
        *fLog << err << "MCerPhotEvt not found... aborting." << endl;
        return kFALSE;
    }

    fCam = (MGeomCam*)pList->FindObject("MGeomCam");
    if (!fCam)
    {
        *fLog << err << "MGeomCam not found... aborting." << endl;
        return kFALSE;
    }

    fTme = 0;
    if (fTimeWindow>0)
    {
        fTme = (MArrivalTime*)pList->FindObject("MArrivalTime");
        if (!fTme)
        {
            *fLog << err << "MArrivalTime not found... aborting." << endl;
            return kFALSE;
        }
    }

    memset(fCut, 0, sizeof(fCut));

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Evaluate software trigger
//
Int_t MFSoftwareTrigger::Process()
{
    switch (fType)
    {
    case kSinglePixelNeighbors:
        fResult = SwTrigger();
        break;
    case kAnyPattern:
        fResult = ClusterTrigger();
        break;
    }

    fCut[fResult ? 0 : 1]++;
    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Prints some statistics about the Basic selections.
//
Int_t MFSoftwareTrigger::PostProcess()
{
    if (GetNumExecutions()==0)
        return kTRUE;

    TString type;
    switch (fType)
    {
    case kSinglePixelNeighbors:
        type = " single pixel trigger";
        break;
    case kAnyPattern:
        type = " any pattern trigger";
        break;
    }

    *fLog << inf << endl;
    *fLog << GetDescriptor() << " execution statistics:" << endl;
    *fLog << " Threshold=" << fThreshold << ", Number=" << (int)fNumNeighbors;
    if (fTimeWindow>0)
        *fLog << ", Time Window=" << fTimeWindow;
    *fLog << endl;
    *fLog << dec << setfill(' ');

    *fLog << " " << setw(7) << fCut[0] << " (" << setw(3) ;
    *fLog << (int)(fCut[0]*100/GetNumExecutions());
    *fLog << "%) Evts fullfilled" << type << endl;
    *fLog << " " << setw(7) << fCut[1] << " (" << setw(3) ;
    *fLog << (int)(fCut[1]*100/GetNumExecutions());
    *fLog << "%) Evts didn't fullfill" << type << "."  << endl;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Read the setup from a TEnv, eg:
//   MFSoftwareTrigger.Threshold:    5
//   MFSoftwareTrigger.NumNeighbors: 4
//   MFSoftwareTrigger.TimeWindow: 3.3
//   MFSoftwareTrigger.TriggerType: SinglePixel, AnyPattern
//
// To switch off time-coincidence use
//   MFSoftwareTrigger.TimeWindow: -1
//
Int_t MFSoftwareTrigger::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    Bool_t rc = kFALSE;
    if (IsEnvDefined(env, prefix, "Threshold", print))
    {
        rc = kTRUE;
        SetThreshold(GetEnvValue(env, prefix, "Threshold", fThreshold));
    }
    if (IsEnvDefined(env, prefix, "NumNeighbors", print))
    {
        rc = kTRUE;
        SetNumNeighbors(GetEnvValue(env, prefix, "NumNeighbors", fNumNeighbors));
    }
    if (IsEnvDefined(env, prefix, "TimeWindow", print))
    {
        rc = kTRUE;
        SetTimeWindow(GetEnvValue(env, prefix, "TimeWindow", fTimeWindow));
    }

    if (IsEnvDefined(env, prefix, "TriggerType", print))
    {
        TString dat = GetEnvValue(env, prefix, "TriggerType", "");
        dat = dat.Strip(TString::kBoth);
        dat.ToLower();

        if (dat == (TString)"singlepixel")
            SetTriggerType(kSinglePixelNeighbors);
        if (dat == (TString)"anypattern")
            SetTriggerType(kAnyPattern);

        rc = kTRUE;
    }

    return rc;
}