////////////////////////////////////////////////////////////////////////////
// Input Containers:
//   MGeomCam
//   MCerPhotEvt
//   MPedestalCam
//
//  Output Containers:
//   MIslands
//
/////////////////////////////////////////////////////////////////////////////
#include "MIslandCalc.h"

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

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

#include "MIslands.h"

#include "MParList.h"

//#include "MImgCleanStd.h"

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

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

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

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


ClassImp(MIslandCalc);


using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor.
//
MIslandCalc::MIslandCalc(const char *name, const char *title)    
  : fIsl(NULL)
{
    fName  = name  ? name  : "MIslandCalc";
    fTitle = title ? title : "Calculate island parameters";
}



Int_t MIslandCalc::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;
      }
    
    
    fIsl = (MIslands*)pList->FindCreateObj(AddSerialNumber("MIslands"));
    if (!fIsl)
      return kFALSE;

    return kTRUE;
}


//tries to clean a bit more the image avoiding noisy islands
 
Int_t MIslandCalc::Process()
{

  Float_t  noise;
  Float_t  signal;

  Int_t    npix = 577;
  Int_t    sflag;
  Int_t    control;
  
  Int_t    nvect = 0;
  Int_t    fIslNum = 0;
  Float_t  fSigToNoise[30];

  Int_t    vect[100][577];
  Int_t    kk[577];
  
  for(Int_t m = 0; m < 50 ; m++)
    for(Int_t n = 0; n < 577 ; n++)
      vect[m][n] = 0;
  
  for(Int_t n = 0; n < 577 ; n++)
    kk[n] = 0;
   
  // loop over all pixels
  for(Int_t idx=0 ; idx<npix ; idx++)
    { 
      const MGeomPix &gpix  = (*fCam)[idx];
      const Int_t    nnmax  = gpix.GetNumNeighbors();
      
      if(idx<0 || idx==npix)
	{
	  cout << "Pixel (software) index out of range [0-576]. Available number of pixels=" << npix << endl;
	  return 1;
	}
      
   
      if(fEvt->IsPixelCore(idx) || fEvt->IsPixelUsed(idx)) 
     	{
	  sflag = 0;
	  
	  for(Int_t j=0; j < nnmax ; j++)
	    {
	      const Int_t idx2 = gpix.GetNeighbor(j);
	      
	      if (idx2 < idx)
		{
		  for(Int_t k = 1; k <= nvect; k++)
		    {
		      if (vect[k][idx2] == 1)
			{
			  sflag = 1;
			  vect[k][idx] = 1;
			}
		    }
		}
	    }
	  
	  if (sflag == 0)
	    {
	      nvect++;
	      vect[nvect][idx] = 1;	     
	    }
	  
	}
    }
  
  fIslNum = nvect;
  
  // Repeated Chain Corrections
  for(Int_t i = 1; i <= nvect; i++)
    {
      for(Int_t j = i+1; j <= nvect; j++)
	{
	  control = 0;
	  for(Int_t k = 0; k < npix; k++)
	    {
	      if (vect[i][k] == 1 && vect[j][k] == 1)
		{
		  control = 1; 
		  //cout << " two vectors coincide... "<< endl;
		  break;
		}
	    }
	  if (control == 1)
	    {
	      for(Int_t k = 0; k < npix; k++)
		{
		  if(vect[j][k] == 1)
		    vect[i][k] = 1;
		  vect[j][k] = 0;
		}	
	      fIslNum = fIslNum-1;	    
	    }
	}
    }
  
  cout << " Let me see..." << endl;
  cout << " Number of vectors: " << nvect << endl;
  cout << " Number of islands: " << fIslNum << endl;      
  
  //examine each island...

  Int_t fPixNum[30];

  for(Int_t i = 1; i<=fIslNum ; i++)
    {
      signal = 0;
      noise = 0;
      fPixNum[i] = 0;
      for(Int_t idx=0 ; idx<npix ; idx++)
	{
	  const MCerPhotPix &pix = (*fEvt)[idx];
	  const MPedestalPix &ped  = (*fPed)[idx];

	  if (vect[i][idx]==1){

	    fPixNum[i]++;
	    signal += pix.GetNumPhotons() * (fCam->GetPixRatio(idx));
	    noise += pow(ped.GetPedestalRms(),2);

	  }	
	}

      cout << "Signal: " << signal<<endl;
      cout << "Noise: " << noise << endl;  
      fSigToNoise[i] = (Float_t)signal/(Float_t)sqrt(noise);

      fIsl->SetIslNum(fIslNum);
      fIsl->SetPixNum(i,fPixNum[i]);
      fIsl->SetSigToNoise(i,fSigToNoise[i]);

      cout << " Island Num " << i << endl;
      cout << " Pixel Number in the island: " << fPixNum[i] << endl;
      cout << " Signal to Noise of the island : " << fSigToNoise[i] << endl;
    	
    }

  return 1;
}