/* ======================================================================== *\
!
!
!   Author(s): Ester Aliu, 2/2004 <aliu@ifae.es>
!  
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MIslandCalc
//
// The Island Calc task calculates the some islands parameters for each 
// of the events such as:
// 
//   - fIslNum             //  number of islands
//   - fIslId[577]         //  island Id
//   - fPixNum[20];        //  number of pixels in the island
//   - fSigToNoise[20];    //  signal to noise of the island
//   - fTimeHi[20][577]    //  mean of the arrival time in the hi gain  
//   - fTimeLo[20][577]    //  mean of the arrival time in the lo gain
//   - fTimeSpreadHi[20];  //  mean arrival time spread of the island  
//   - fTimeSpreadLo[20];  //  mean arrival time spread of the island
//
//
// Input Containers:
//   MGeomCam
//   MCerPhotEvt
//   MPedestalCam
//   MArrivalTimeCam
//
//  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 "MGeomPix.h"
#include "MGeomCam.h"

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

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

#include "MArrivalTimeCam.h"
#include "MArrivalTimePix.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;
      }

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


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;
 
  Int_t    vect[50][577];

  Int_t    zeros[50];
  
  for(Int_t m = 0; m < 50 ; m++)
    for(Int_t n = 0; n < npix ; n++)
	vect[m][n] = 0;
    
  for(Int_t n = 0; n < 50 ; n++)
    zeros[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->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; 
		  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;
		  zeros[j] = 1;
		}	
	      fIslNum = fIslNum-1;	    
	    }
	  
	}
    }
  
 
  Int_t l = 1;

  for(Int_t i = 1;  i<= nvect ; i++)
    {
      if (zeros[i] == 0)
	{
	  for(Int_t k = 0; k<npix; k++)
	    {
	      vect[l][k] = vect[i][k];
	    }
	  l++;
	}
    }
   
  
  //set the number of islands in one event
  fIsl->SetIslNum(fIslNum);

  //examine each island...
  Int_t fPixNum[fIslNum];
  Float_t fSigToNoise[fIslNum];
  Float_t time[577];
  Float_t timeVariance[fIslNum];
  //Float_t timeHi[577], timeLo[577];  
  //Float_t timeVarianceHi[fIslNum];
  //Float_t timeVarianceLo[fIslNum];
 
  //reset the "sets" functions
  if (fIslNum <1)
    fIsl->SetIslNum(0);

  for(Int_t i = 0; i<20 ;i++)
    {
      fIsl->SetPixNum(i,-1);
      fIsl->SetSigToNoise(i,-1);
      fIsl->SetTimeSpread(i,-1);
      //      fIsl->SetTimeSpreadHi(i,-1);
      // fIsl->SetTimeSpreadLo(i,-1);

      for(Int_t idx = 0; idx<npix; idx++)
	{
	  fIsl->SetIslId(idx, -1);
	  fIsl->SetArrivalTime(i, idx, -1 );
	  // fIsl->SetArrivalTimeHiGain(i, idx, -1 );
	  // fIsl->SetArrivalTimeLoGain(i, idx, -1);
	}
    }
 
    
  for(Int_t i = 1; i<=fIslNum ; i++)
    {
      Int_t n = 0;
      Int_t ncore = 0;

      Float_t MINhi = 10000;
      Float_t MINlo = 10000;
      Float_t MAXhi = 0;
      Float_t MAXlo = 0;
      Float_t MIN = 10000;
      Float_t MAX = 0;

      //cout << "Isl #" << i << endl;

      signal = 0;
      noise = 0;
      fPixNum[i-1] = 0;
      timeVariance[i-1] = 0;
      // timeVarianceHi[i-1] = 0;
      // timeVarianceLo[i-1] = 0;

      for(Int_t idx=0 ; idx<npix ; idx++)
	{
	  const MCerPhotPix &pix = (*fEvt)[idx];
	  const MPedestalPix &ped  = (*fPed)[idx];
	  const MArrivalTimePix &timepix = (*fTime)[idx];
	  
	  if (vect[i][idx]==1){
    
	    fPixNum[i-1]++;
	    signal += pix.GetNumPhotons() * (fCam->GetPixRatio(idx));
	    noise += pow(ped.GetPedestalRms(),2);

	    time[n] = timepix.IsLoGainUsed() ? timepix.GetArrivalTimeLoGain() : timepix.GetArrivalTimeHiGain();

	    // timeHi[n] = timepix.GetArrivalTimeHiGain(); 
	    // cout << "timeHi" << timeHi[n] << endl;	    

	    //timeLo[n] = timepix.GetArrivalTimeLoGain(); 
	    //cout << "timeHi" << timeLo[n] << endl;	    
	
	    if (fEvt->IsPixelCore(idx)){ 
	      
	      /*  if (timeHi[n] < MINhi)
		MINhi = timeHi[n];
	      if (timeLo[n] < MINlo)
		MINlo = timeLo[n];

	      if (timeHi[n] > MAXhi)
		MAXhi = timeHi[n];
	      if (timeLo[n] > MAXlo)
		MAXlo = timeLo[n];

	      */

	      if (time[n] > MAX)
		MAX = time[n];
	      if (time[n] < MIN)
		MIN = time[n];
	      
	      //control2[n] = 1;   
	      //timeVarianceHi[i-1] += timeHi[n];
	      //timeVarianceLo[i-1] += timeLo[n];
	      ncore++;
	    }

	    fIsl->SetIslId(idx, i-1);
	    fIsl->SetArrivalTime(i-1, idx, time[n]);
	    // fIsl->SetArrivalTimeHiGain(i-1, idx, timeHi[n]);
	    // fIsl->SetArrivalTimeLoGain(i-1, idx, timeLo[n]);
	
	    n++;
	  }
	  
	}

      Float_t mean = timeVariance[i-1]/ncore;
      // Float_t meanHi = timeVarianceHi[i-1]/ncore;
      // Float_t meanLo = timeVarianceLo[i-1]/ncore;
      
      timeVariance[i-1] = 0;
      //timeVarianceHi[i-1] = 0;
      //timeVarianceLo[i-1] = 0;

      /*
      //old
	for (Int_t k = 0; k <n ; k++)
	{ 
	if (control2[k] == 1){
	
	timeVarianceHi[i-1] += pow(timeHi[k]- meanHi,2); 
	timeVarianceLo[i-1] += pow(timeLo[k]- meanLo,2); 
	}
	}
	timeVarianceHi[i-1] = sqrt(timeVarianceHi[i-1]/(ncore-1)); 
	timeVarianceLo[i-1] = sqrt(timeVarianceLo[i-1]/(ncore-1)); */
      
      //timeVarianceHi[i-1] = (MAXhi - MINhi)/ncore; 
      //timeVarianceLo[i-1] = (MAXlo - MINlo)/ncore; 

      timeVariance[i-1] = (MAX - MIN)/ncore; 
      timeVariance[i-1] = (MAX - MIN)/ncore; 

      fSigToNoise[i-1]= (Float_t)signal/(Float_t)sqrt(noise);

      fIsl->SetPixNum(i-1,fPixNum[i-1]);
      fIsl->SetSigToNoise(i-1,fSigToNoise[i-1]);
      fIsl->SetTimeSpread(i-1, timeVariance[i-1]);
      //fIsl->SetTimeSpreadHi(i-1, timeVarianceHi[i-1]);
      //fIsl->SetTimeSpreadLo(i-1, timeVarianceLo[i-1]);

      //cout << " TimeHi Spread: " << timeVarianceHi[i-1] << endl;
      //cout << " # core pixels: " << ncore << endl;
      //cout << " # used pixels: " << n << endl;
      //cout << " SigToNoise: " << fSigToNoise[i-1] << endl << endl;
      
    } 
  
  fIsl->SetReadyToSave();
  
  return 1;
  
}