source: trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsClean.cc@ 5295

/* ======================================================================== *\
!
!
!   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 "MIslandsClean.h"

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

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

#include "MIslands.h"
#include "MImgIsland.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(MIslandsClean);


using namespace std;

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


Int_t MIslandsClean::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 MIslandsClean::Process()
{
  //
  //eliminates the island with a signal-to-noise 
  //lower than a given limit 
  //
  //if ( fIslandCleaningMethod == kNoTiming ){
  
  MImgIsland *imgIsl = new MImgIsland;
  TIter Next(fIsl->GetList());

  Int_t pixNum = 0;  
  Int_t idPix = -1;

  if ( fIslandCleaningMethod == 1 ) {
    while ((imgIsl=(MImgIsland*)Next())) {
      
      if(imgIsl->GetSigToNoise() < fIslCleanThres)
        {
          pixNum = imgIsl->GetPixNum();
          
          for(Int_t k = 0; k<pixNum; k++)
            {
              idPix = imgIsl->GetPixList(k);
              MCerPhotPix &pix  = (*fEvt)[idPix];
              pix.SetPixelUnused();
            }
        }
    }   
  }  
  
  //
  //eliminates the island with a time spread  
  //higher than a given limit 
  //
  //else if( fIslandCleaningMethod == kTiming ){
  else if( fIslandCleaningMethod == 0 ){
    while ((imgIsl=(MImgIsland*)Next())) {

      //fIslCleanThreshold has different values, FIXME, put two variables
      if(imgIsl->GetTimeSpread() > fIslCleanThres)
        {
          pixNum = imgIsl->GetPixNum();
          
          for(Int_t k = 0; k<pixNum; k++)
            {
              idPix = imgIsl->GetPixList(k);
              MCerPhotPix &pix  = (*fEvt)[idPix];
              pix.SetPixelUnused();
            }
        }
    }
    
  }
  
  //
  //eliminates all the islands except the continent, 
  //i.e. the larger island in the event
  //
  else if( fIslandCleaningMethod == 3 ){
    Int_t i = 0;
    while ((imgIsl=(MImgIsland*)Next())) {
      if (i != 0){
        
        pixNum = imgIsl->GetPixNum();
        
        for(Int_t k = 0; k<pixNum; k++)
          {
            idPix = imgIsl->GetPixList(k);
            MCerPhotPix &pix  = (*fEvt)[idPix];
            pix.SetPixelUnused();
          }
      }
      i++;   
    }
  }
  
  fEvt->SetReadyToSave();
  
  return kTRUE;
  
}