/////////////////////////////////////////////////////////////////////////////
//
// EEezCluster
// 
// Author: Jason C. Webb <jwebb@iucf.indiana.edu>
//
// Storage class for clusters.  Inherits from EEezTower.  Will have a 
// method to build clusters from a specified seed tower.  Is meant to
// be a persistent object... i.e. an array of 720 m_Clusters sitting
// within the cluster maker.  Neighboring clusters can be easily associated
// with this cluster by storing pointers (much like was done for the towers).
//
/////////////////////////////////////////////////////////////////////////////

#include "EEezCluster.h"

#include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
#include "StEEmcUtil/StEEmcSmd/EEmcSmdGeom.h"

#define DEBUG 0

/////////////////////////////////////////////////////////////////////////////
ClassImp(EEezCluster);

 EEezCluster::EEezCluster() 
{
  // Class initialization.
  
  m_EETowGeom = new EEmcGeomSimple();
  m_EESmdGeom = EEmcSmdGeom::instance();
  m_EnergyCut = 0.;

  for ( Int_t i = 0; i < 2; i++ ) {
    m_UFirst[i] = 0;
    m_VFirst[i] = 0;
    m_ULast[i]  = 0;
    m_VLast[i]  = 0;
  }

}

/////////////////////////////////////////////////////////////////////////////

 Int_t EEezCluster::hasTower ( EEezTower *tower ) 
{
  // Returns true if the tower has been associated with this cluster
  // NOTE: should add a pointer to EEezTower to hold all clusters 
  // which it has been associated with.
  

  Int_t dEta = tower -> getEtabin() - m_SeedTower -> getEtabin();
  Int_t dPhi = tower -> getPhibin() - m_SeedTower -> getPhibin();

  Int_t ret =
    ( (dEta*dEta)<=1 ) &&
    ( (dPhi*dPhi)<=1 ) &&
    tower -> getEnergy() > m_EnergyCut;

  return ret;

}

/////////////////////////////////////////////////////////////////////////////

 void EEezCluster::buildCluster( EEezTower *seed ) 
{ 
  // Performs clustering of the specified seed tower.  The energies
  // and momentum (energy times direction from vertex) of the seed's 
  // neighbors are summed, if the neighboring energy exceeds a
  // specified threshold.
  

  // First, set the seed tower for this cluster
  m_SeedTower = seed;
  m_NumTowers = 1;
  m_NumSeeds  = 1;
  m_Index     = seed -> getIndex();

  //-- Add this seed tower to the list of seed towers 
  //-- for this cluster (generalized cluster may have
  //-- more than one seed tower).
  m_SeedTowers.push_back(seed);


  // Get the list of neighboring towers, and the list of
  // neighboring seed towers
  EEezTowerPtrVec_t fNeighbors     = seed -> getNeighbors();
  EEezTowerPtrVec_t fNeighborSeeds = seed -> getNeighborSeeds();
  m_NumSeeds += fNeighborSeeds.size();


  // Extract seed-tower parameters to local varaibles  
//   Int_t   fIndex  = seed -> getIndex();  // Tower index (0-719)
//   Int_t   fEta    = fIndex % 12;         // Etabin 
//   Int_t   fPhi    = fIndex / 12;         // Phibin
//   Int_t   fSec    = fPhi / 5;            // Sector
//   Int_t   fSub    = fPhi % 5;            // Subsector

  Int_t fIndex = seed -> getIndex();  // Tower index (0-719)
  Int_t fEta   = seed -> getEtabin(); // Etabin  (0-11)
  Int_t fPhi   = seed -> getPhibin(); // Phibin  (0-59)
  Int_t fSec   = seed -> getSector(); // Sector  (0-11)
  Int_t fSub   = seed -> getSubSector(); // SubSector (0-4)
  Float_t fEnergy = seed -> getEnergy(); // Energy of seed tower
  
  // Intial momentum vector
  TVector3 fMomentum = fEnergy * ( m_EETowGeom -> getTowerCenter(fSec,fSub,fEta).Unit() );


  //
  // Initialize pointers to first and last U and V strips
  //
  m_UFirst[0] = seed -> getUStrip(0);
  m_VFirst[0] = seed -> getVStrip(0);
  Int_t nu = seed -> getNUStrips();
  Int_t nv = seed -> getNVStrips();
  m_ULast[0] = seed -> getUStrip(nu-1);
  m_VLast[0] = seed -> getVStrip(nv-1);


  //
  // Our initial cluster model will be the seed tower plus the 
  // 8 surrounding towers.  Futher analysis may bring the SMD
  // in to refine the cluster.  For overlapping clusters, the
  // energy is split between the two, proportional to the 
  // energies of each seed.
  //


  // Iterate over neighboring towers and calculate energy,
  // eta and phi.
  EEezTowerPtrVecIter_t iter = fNeighbors.begin();

  // Running sum for energy.  Eta and phi will be computed from
  // the resulting summed momentum vector.
  Float_t fSumEnergy = fEnergy;
  

  // Loop over all neighboring towers
  while ( iter != fNeighbors.end() ) { 
    
    // Get the index of this tower, etc...
    Float_t fEnergy1 = (*iter)->getEnergy(); // Energy of this tower
    if ( fEnergy1 < m_EnergyCut ) continue;
    Int_t fIndex1  = (*iter)->getIndex();  // Tower index
//     Int_t fEta1    = fIndex1 % 12;         // Etabin
//     Int_t fPhi1    = fIndex1 / 12;         // Phibin
//     Int_t fSec1    = fPhi1 / 5;            // Sector
//     Int_t fSub1    = fPhi1 % 5;            // Subsector 
    Int_t fEta1 = (*iter)->getEtabin();
    Int_t fPhi1 = (*iter)->getPhibin();
    Int_t fSec1 = (*iter)->getSector();
    Int_t fSub1 = (*iter)->getSubSector(); 

    // Increment our tower count
    m_NumTowers++;

    // Get the momentum vector pointing at this tower
    TVector3 fMomentum1 = fEnergy1 * ( m_EETowGeom -> getTowerCenter(fSec1,fSub1,fEta1).Unit() );

    // Get the energy summed over all adjacent seed towers
    // for the current neighboring tower.  The weight for 
    // this tower is equal to the energy of the seed tower
    // divided by this number.

    EEezTowerPtrVec_t     fAdjSeeds = (*iter) -> getNeighborSeeds();
    EEezTowerPtrVecIter_t iterAdjSeeds = fAdjSeeds.begin();
    Float_t               fEnergyAdjSeeds = 0.;
    while ( iterAdjSeeds != fAdjSeeds.end() ) { 
      fEnergyAdjSeeds += (*iterAdjSeeds)->getEnergy(0);
      iterAdjSeeds++;
    }
    Float_t fWeight = seed -> getEnergy();
    if ( fEnergyAdjSeeds > 0. ) {
      fWeight /= fEnergyAdjSeeds;
    }
    else {
      fWeight = 1.0;
    }

    // Do the math
    fSumEnergy += ( fEnergy1   * fWeight );
    fMomentum  += ( fMomentum1 * fWeight );


    //--
    //-- Get pointers to the first and last U and V Strips
    //-- contained w/in the fiducial volume of this tower
    //--
    nu = (*iter) -> getNUStrips();
    nv = (*iter) -> getNVStrips();
    EEezStrip *uf = (*iter) -> getUStrip(0);
    EEezStrip *vf = (*iter) -> getVStrip(0);
    EEezStrip *ul = (*iter) -> getUStrip(nu-1);
    EEezStrip *vl = (*iter) -> getVStrip(nv-1);

    //-- Towers in the same sector with seed
    if ( fSec == fSec1 ) {
      if ( m_UFirst[0] -> getIndex() > uf -> getIndex() )
	m_UFirst[0] = uf;
      if ( m_VFirst[0] -> getIndex() > vf -> getIndex() )
	m_VFirst[0] = vf;
      if ( m_ULast[0] -> getIndex() < ul -> getIndex() )
	m_ULast[0] = ul;
      if ( m_VLast[0] -> getIndex() < vl -> getIndex() )
	m_VLast[0] = vl;      
    }
    //-- Towers in sector other than seed
    else {
      if ( !m_UFirst[1] ) {
	m_UFirst[1] = uf;
	m_VFirst[1] = vf;
	m_ULast[1]  = ul;
	m_VLast[1]  = vl;
      }
      if ( m_UFirst[1] -> getIndex() > uf -> getIndex() )
	m_UFirst[1] = uf;
      if ( m_VFirst[1] -> getIndex() > vf -> getIndex() )
	m_VFirst[1] = vf;
      if ( m_ULast[1] -> getIndex() < ul -> getIndex() )
	m_ULast[1] = ul;
      if ( m_VLast[1] -> getIndex() < vl -> getIndex() )
	m_VLast[1] = vl;    
    }



    //-- Next neighboring tower
    iter++;

  }



  // NOTE TO SELF-- we are NOT doing a weighted average here!  
  m_ClusterEnergy = fSumEnergy;



#if DEBUG
  std::cout << " Final: " 
	    << fMomentum[0] << " "
 	    << fMomentum[1] << " "
 	    << fMomentum[2] << " "
 	    << std::endl;

  std::cout << " Energy: " 
 	    << m_ClusterEnergy << " "
 	    << fMomentum.Mag() << std::endl;

  std::cout << std::endl;
#endif

  

  // There is a ~0.1% difference between the weighted energy sum,
  // and the magnitude of the resulting momentum vector.  Rescale
  // the momentum vector by our calculated energy.
  
  m_ClusterMomentum = m_ClusterEnergy * ( fMomentum.Unit() );


}


/////////////////////////////////////////////////////////////////////////////

 void EEezCluster::setSeedTower( EEezTower *tower ) 
{
  // Sets the seed tower of this cluster... I don't see the need for
  // this method any longer... will probably delete it.
  
  
  m_SeedTower = tower;
  
}

/////////////////////////////////////////////////////////////////////////////

 void EEezCluster::Consume ( EEezCluster *cluster ) 
{ 
  // The idea behind this method was to provide a way to merge clusters
  // when doing mixed-event backgrounds.  At the momentum, this method
  // does very little and should not be used.
  

  // Simple for now...
 
  m_ClusterMomentum += cluster -> getMomentum();
  m_ClusterEnergy   += cluster -> getEnergy();


}

/////////////////////////////////////////////////////////////////////////////

 void EEezCluster::setType( EEezClusterType type )
{
  // Sets the type of the cluster.
  //
  // EEezCluster::kEEezClusterSingle
  //
  // "Single Seed" cluster.  The cluster is comprised of 8 towers
  // surrounding a single seed tower, such that E_seed / E_cluster
  // is larger than the value specified by setShapeLimit().
  //
  // EEezCluster::kEEezClusterDouble (not implemented)
  //
  // "Double Seed" cluster.  Single-seed clusters which fail the
  // shape limit may be classified as a double-seed cluster.  The
  // cluster will be comprised of the 10 neighbors of two non-
  // diagonally adjacent seed towers.  It is required that the sum
  // E_seed1 + E_seed2 / E_cluster is larger than the value 
  // specified by setShapeLimit().
  //
  // EEezCluster::kEEezClusterDiagonal (not implemented)
  //
  // "Double Seed Diagonal" cluster.  The cluster is comprised of
  // the 12 towers surrounding two diagonally-adjacent seed towers.
  // The sum of the two seeds is required to exceed the value
  // specified by setShapeLimit().
  //
  // EEezCluster::kEEezClusterOther (not implemented)
  //
  // This is a catch-all category which will be used to accumulate
  // multiple-adjacent seed towers, such as may be frequently encountered
  // in AuAu running.   These will be stored in a seperate vector
  // from "normal" clusters.... This type will be better defined
  // once implemented.
  //
  mType = type;

}
 Int_t EEezCluster::getType() 
{
  // Returns the type of cluster.  See setType() for definitions.
  return mType;
}