source: trunk/Mars/mimage/MImgCleanFAMOUS.cc@ 20003

/* ======================================================================== *\
!
! *
! * 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): Thomas Bretz, 06/2020 <mailto:tbretz@physik.rwth-aachen.de>
!
!   Copyright: MAGIC Software Development, 2000-2020
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MImgCleanFAMOUS
//
// This image cleaning algorithm was introduced by J. Audehm (Master
// thesis, RWTH Aachen, 2020) based on a cleaning by M. Schaufel
// (Master thesus, RWTH Aachen, 2017). It was slightly modified
// to fit into the framework.
//
/////////////////////////////////////////////////////////////////////////////
#include "MImgCleanFAMOUS.h"

#include <TEnv.h>

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

#include "MParList.h"

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

#include "MSignalPix.h"
#include "MSignalCam.h"

ClassImp(MImgCleanFAMOUS);

using namespace std;

static const TString gsDefName  = "MImgCleanFAMOUS";
static const TString gsDefTitle = "Task to perform image cleaning";

const TString MImgCleanFAMOUS::gsNameSignalCam ="MSignalCam";  // default name of the 'MSignalCam' container
const TString MImgCleanFAMOUS::gsNameGeomCam   ="MGeomCam";    // default name of the 'MGeomCam' container

// --------------------------------------------------------------------------
//
// Default constructor. Here you can specify the cleaning method and levels. 
// If you don't specify them the 'common standard' values 3.0 and 2.5 (sigma
// above mean) are used.
// Here you can also specify how many rings around the core pixels you want 
// to analyze (with the fixed lvl2). The default value for "rings" is 1.
//
MImgCleanFAMOUS::MImgCleanFAMOUS(const char *name, const char *title)
    : fSignalMin(40), fSignalMax(4000), fTimeMin(110), fTimeMax(150),
    fSlopeMin(0), fSlopeMax(100), fTimeWindow(5),
    fNameGeomCam(gsNameGeomCam), fNameSignalCam(gsNameSignalCam)
{
    fName  = name  ? name  : gsDefName.Data();
    fTitle = title ? title : gsDefTitle.Data();
}

void MImgCleanFAMOUS::ResetCleaning() const
{
    //
    // check the number of all pixels against the noise level and
    // set them to 'unused' state if necessary
    //
    const UInt_t npixevt = fEvt->GetNumPixels();
    for (UInt_t idx=0; idx<npixevt; idx++)
    {
        MSignalPix &pix = (*fEvt)[idx];
        if (pix.IsPixelUnmapped())
            continue;

        pix.SetPixelUnused();
        pix.SetPixelCore(kFALSE);
        pix.SetIdxIsland(-1);
    }
}

// --------------------------------------------------------------------------
//
void MImgCleanFAMOUS::FindIsland(Int_t idx) const
{
    // Get the pixel information of a pixel with this index
    MSignalPix &pix = (*fEvt)[idx];
    if (!pix.IsPixelUsed())
        return;

    // Get geometrical description of pixel idx
    const MGeom &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++)
    {
        // Index of neighbor pixel
        const UInt_t ipix = gpix.GetNeighbor(i);

        // Signal of neighbor pixel
        const MSignalPix &npix = (*fEvt)[ipix];

        if (TMath::Abs(npix.GetArrivalTime()-pix.GetArrivalTime())>=fTimeWindow)
            continue;

        if (npix.GetNumPhotons()<=0)
            continue;

        // FIXME: Checks on minimum signal and slope parameter?

        if (npix.IsPixelUnmapped())
            continue;

        pix.SetPixelUsed();
        pix.SetIdxIsland(1);

        FindIsland(ipix);
    }
}


void MImgCleanFAMOUS::DoCleaning() const
{
    UInt_t max_idx = -1;

    // Serach for the brightest none-saturating pixel
    // that is consistent with originating form a shower
    const UInt_t npixevt = fEvt->GetNumPixels();
    for (UInt_t idx=0; idx<npixevt; idx++)
    {
        MSignalPix &pix = (*fEvt)[idx];

        // Check that the slopd of the rising edge looks reasonable
        if (pix.GetTimeSlope() <= fSlopeMin || pix.GetTimeSlope()>fSlopeMax)
            continue;

        // Check that pixel aligns with the trigger position
        if (pix.GetArrivalTime() <= fTimeMin || pix.GetArrivalTime()>fTimeMax)
            continue;

        // Check if this is above the cleaning level and not saturating
        if (pix.GetNumPhotons() <= fSignalMin || pix.GetNumPhotons()>fSignalMax)
            continue;

        // Ignore unmapped pixels
        if (pix.IsPixelUnmapped())
            continue;

        pix.SetPixelCore();

        if (max_idx<0 || pix.GetNumPhotons()>(*fEvt)[max_idx].GetNumPhotons())
            max_idx = idx;
    }

    if (max_idx<0)
        return;

    (*fEvt)[max_idx].SetPixelUsed();
    FindIsland(max_idx);
}

// --------------------------------------------------------------------------
//
//  Check if MEvtHeader exists in the Parameter list already.
//  if not create one and add them to the list
//
Int_t MImgCleanFAMOUS::PreProcess(MParList *pList)
{
    fCam = (MGeomCam*)pList->FindObject(AddSerialNumber(fNameGeomCam), "MGeomCam");
    if (!fCam)
    {
        *fLog << err << fNameGeomCam << " [MGeomCam] not found (no geometry information available)... aborting." << endl;
        return kFALSE;
    }
    fEvt = (MSignalCam*)pList->FindObject(AddSerialNumber(fNameSignalCam), "MSignalCam");
    if (!fEvt)
    {
        *fLog << err << fNameSignalCam << " [MSignalCam] not found... aborting." << endl;
        return kFALSE;
    }

    Print();

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Cleans the image.
//
Int_t MImgCleanFAMOUS::Process()
{
    ResetCleaning();
    DoCleaning();

    // Takes roughly 10% of the time
    fEvt->CalcIslands(*fCam);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Print descriptor and cleaning levels.
//
void MImgCleanFAMOUS::Print(Option_t *o) const
{
    *fLog << all << GetDescriptor() << " using" << endl;
    *fLog << " * " << fSignalMin << " < Signal <= " << fSignalMax << endl;
    *fLog << " * " << fTimeMin   << " < Time <= "   << fTimeMax   << endl;
    *fLog << " * " << fSlopeMin  << " < Slope <= "  << fSlopeMax  << endl;
    *fLog << " * " << "Delta T < "  << fTimeWindow << endl;
}

// --------------------------------------------------------------------------
//
// Read the setup from a TEnv, eg:
//   MImgCleanFAMOUS.SignalMin: 40
//   MImgCleanFAMOUS.SignalMax: 4000
//   MImgCleanFAMOUS.TimeMin: 110
//   MImgCleanFAMOUS.TimeMax: 150
//   MImgCleanFAMOUS.SlopeMin: 0
//   MImgCleanFAMOUS.SlopeMax: 100
//   MImgCleanFAMOUS.TimeWindow: 5
//
Int_t MImgCleanFAMOUS::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    Bool_t rc = kFALSE;
    if (IsEnvDefined(env, prefix, "SignalMin", print))
    {
        rc = kTRUE;
        fSignalMin = GetEnvValue(env, prefix, "SignalMin", fSignalMin);
    }
    if (IsEnvDefined(env, prefix, "SignalMax", print))
    {
        rc = kTRUE;
        fSignalMax = GetEnvValue(env, prefix, "SignalMax", fSignalMax);
    }
    if (IsEnvDefined(env, prefix, "TimeMin", print))
    {
        rc = kTRUE;
        fTimeMin = GetEnvValue(env, prefix, "TimeMin", fTimeMin);
    }
    if (IsEnvDefined(env, prefix, "TimeMax", print))
    {
        rc = kTRUE;
        fTimeMax = GetEnvValue(env, prefix, "TimeMax", fTimeMax);
    }
    if (IsEnvDefined(env, prefix, "SlopeMin", print))
    {
        rc = kTRUE;
        fSlopeMin = GetEnvValue(env, prefix, "SlopeMin", fSlopeMin);
    }
    if (IsEnvDefined(env, prefix, "SlopeMax", print))
    {
        rc = kTRUE;
        fSlopeMax = GetEnvValue(env, prefix, "SlopeMax", fSlopeMax);
    }
    if (IsEnvDefined(env, prefix, "TimeWindow", print))
    {
        rc = kTRUE;
        fTimeWindow = GetEnvValue(env, prefix, "TimeWindow", fTimeWindow);
    }

    return rc;
}