source: trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandClean.cc@ 4450

/* ======================================================================== *\
!
!
!   Author(s): Ester Aliu, 3/2004
!  
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MIslandClean
//
// The Island Cleaning task selects the islands you use for the Hillas
// parameters calculation, after the normal image cleaning.
//
// There are two methods to make the selection: 
//
//    - No time method, as used is Whipple. It calculates an island parameter 
//   called "signal to noise" and adds a new threshold that eliminates some   
//   of the islands. The way the island is eliminated is seeting the pixels 
//   of the islands as UNUSED 
//
//    - Time method, taking profit of the time information in MAGIC.
//   Calculates the maximum time difference (in time slices) for each island 
//   and corrects for the island size. With an other new threshold "noise" 
//   islands are supposed to be eliminated.     
//
// Other cleanings that are allowed in this code are:
// 
//    - Resting only with the continent, i.e. the larger island
//    
//  Example:
//
//  MIslands   isl;
//  isl.SetName("MIslands1");

//  MImgCleanStd clean;
//  MIslandClean islclean(0.2);
//  islclean.SetInputName("MIslands1");  
//  islclean.SetMethod(0); // for timing method 
//
//  tlist.AddToList(&clean);
//  tlist.AddToList(&islclean);
//
//
//  Input Containers:
//    MGeomCam
//    MCerPhotEvt
//    MPedestalCam
//    MArrivalTime
//    MIslands
//
//  Output Containers:
//    MCerPhotEvt
//
/////////////////////////////////////////////////////////////////////////////
#include "MIslandClean.h"

#include <stdlib.h>       // atof                                         
#include <fstream>        // ofstream, SavePrimitive

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

#include "MIslands.h"

#include "MParList.h"

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

#include "MCerPhotPix.h"
#include "MCerPhotEvt.h"

#include "MPedestalCam.h"
#include "MPedestalPix.h"

#include "MArrivalTimeCam.h"
#include "MArrivalTimePix.h"

ClassImp(MIslandClean);


using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor.
//
MIslandClean::MIslandClean(const Float_t newThres, const char *name, const char *title)    
  : fIsl(NULL), fIslandCleaningMethod(kNoTiming), fIslCleanThres(newThres)
{
    fName  = name  ? name  : "MIslandClean";
    fTitle = title ? title : "Clean islands";
}


Int_t MIslandClean::PreProcess (MParList *pList)
{
    fCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
    if (!fCam)
    {
        *fLog << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl;
        return kFALSE;
    }

    fEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber("MCerPhotEvt"));
    if (!fEvt)
    {
        *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl;
        return kFALSE;
    }

    fPed = (MPedestalCam*)pList->FindObject(AddSerialNumber("MPedestalCam"));
    if (!fPed)
      {
        *fLog << dbginf << "MPedestalCam not found... aborting." << endl;
        return kFALSE;
      }

    fTime = (MArrivalTimeCam*)pList->FindObject(AddSerialNumber("MArrivalTimeCam"));
    if (!fTime)
      {
        *fLog << dbginf << "MArrivalTimeCam not found... aborting." << endl;
        return kFALSE;
      }

    if (strlen(fIslName) > 0)
      fIsl = (MIslands*)pList->FindObject(AddSerialNumber(fIslName));
    else
      fIsl = (MIslands*)pList->FindObject(AddSerialNumber("MIslands"));
    if (!fIsl)
      {
        *fLog << dbginf << "MIslands not found... aborting." << endl;
        return kFALSE;
      }
    
    return kTRUE;
}



Int_t MIslandClean::Process()
{
  //
  //eliminates the island with a signal-to-noise 
  //lower than a given limit 
  //
  //if ( fIslandCleaningMethod == kNoTiming ){
  if ( fIslandCleaningMethod == 1 ){
    Int_t nisl = fIsl->GetIslNum();
    
    for(Int_t isl = 0; isl<nisl ; isl++)
      {
        if(fIsl->GetSigToNoise(isl) < fIslCleanThres)
          {
            for(Int_t idx = 0; idx<577; idx++)
              {
                MCerPhotPix &pix  = (*fEvt)[idx];
                
                if (fIsl->GetIslId(idx) == isl)
                  pix.SetPixelUnused();
              }
          }
      } 
  }  
    
  //
  //eliminates the island with a time spread  
  //higher than a given limit 
  //
  //else if( fIslandCleaningMethod == kTiming ){
  else if( fIslandCleaningMethod == 0 ){
    Int_t nisl = fIsl->GetIslNum();
    
    for(Int_t isl = 0; isl<nisl ; isl++)
      {
        //fIslCleanThreshold has different values, FIXME, put two variables
        
        if(fIsl->GetTimeSpread(isl) > fIslCleanThres)
          {
            for(Int_t idx = 0;idx<577; idx++)
              {
                MCerPhotPix &pix  = (*fEvt)[idx];
                
                if (fIsl->GetIslId(idx) == isl)
                  pix.SetPixelUnused();
              }
          }
        
      }     
  }

  //
  //eliminates all the islands except the continent, 
  //i.e. the larger island in the event
  //
  else if( fIslandCleaningMethod == 3 ){
    Int_t nisl = fIsl->GetIslNum();

    Int_t max = -1000;    
    Int_t idMax;

    for(Int_t isl = 0; isl<nisl ; isl++)
      {
        if (fIsl->GetPixNum(isl)>max)
         {
           max = fIsl->GetPixNum(isl);
           idMax = isl;
         } 
      }     

    for(Int_t isl = 0; isl<nisl ; isl++)
      {
        if (isl != idMax)
          {
            for(Int_t idx = 0;idx<577; idx++)
              {
                MCerPhotPix &pix  = (*fEvt)[idx];
                
                if (fIsl->GetIslId(idx) == isl)
                  pix.SetPixelUnused();
              }
          }
      }
  }
  
  fEvt->SetReadyToSave();
  
  return kTRUE;
  
}