#include "Cluster.h"
#include "McPart.h"
#include "EDeposit.h"
#include "TMath.h"
#include "TObjArray.h"


ClassImp(Cluster)

  //////////////////////////////////////////////////////////////////////////
  //                                                                      
  // Cluster                                                              
  //                                                                      
  // The Cluster class is used to extract cluster information from        
  // the GISMO ascii output file.                                         
  //                                                                      
  //////////////////////////////////////////////////////////////////////////

///_________________________________________________________________________
 Cluster::Cluster():m_energy(0),m_energyTheta(0),m_energyPhi(0),m_energyThetaRMS(0),m_energyPhiRMS(0),m_eR(0),m_eRRMS(0),m_clusterStartDepth(0){
  // Default constructor
}
//_________________________________________________________________________
 Cluster::~Cluster(){
  // Destructor
  Clean();
}
//_________________________________________________________________________
 void Cluster::Clean(){
  // This function cleans up all the TObjArrays
  m_mcEList.Delete();
  m_IndexList.Delete();
  m_mcList.Delete();
  m_hitList.Delete();
}
///_______________________________________________________________________
 Float_t Cluster::GetEnergy(){
  // Returns the total energy in the cluster
Float_t Etot = 0.0;
 if(!m_energy){
   Int_t nentries = MCEList()->GetEntries();
   for(Int_t i=0; i < nentries; i++){
     Etot += ((EDeposit*)MCEList()->At(i))->GetEDeposit();
   }
   return Etot;
 }
 else return m_energy;
}
///_______________________________________________________________________
 Float_t Cluster::GetEnergyTheta(){
  // Returns the energy weighted theta average of the cluster
  Int_t nentries = HitList()->GetEntries();
  Float_t sumEtheta,sumE;
  if(!m_energyTheta){
    for(Int_t i=0; i<nentries; i++){
      Float_t E = ((EDeposit*)MCEList()->At(i))->GetEDeposit();
      Float_t Etheta = ((towerID*)HitList()->At(i))->theta()*E;
      sumEtheta += Etheta;
      sumE += E;
    }
    return(sumEtheta/sumE);
  }
  else return m_energyTheta;
}
///_______________________________________________________________________
 Float_t Cluster::GetEnergyPhi(){
  // Retruns the energy weighted phi average of the cluster
  Int_t nentries = HitList()->GetEntries();
  Float_t sumEphi,sumE;
  if(!m_energyPhi){
    for(Int_t i=0; i<nentries; i++){
      Float_t E = ((EDeposit*)MCEList()->At(i))->GetEDeposit();
      Float_t Ephi = ((towerID*)HitList()->At(i))->phi()*E;
      sumEphi += Ephi;
      sumE += E;
    }
    return(sumEphi/sumE);
  }
  else return m_energyPhi;
}
///_______________________________________________________________________
 Float_t Cluster::GetEnergyR(){
  return m_eR;
}
///_______________________________________________________________________  
 Float_t Cluster::GetEnergyThetaRMS(){
  return m_energyThetaRMS;
}
///_______________________________________________________________________
 Float_t Cluster::GetEnergyPhiRMS(){
  return m_energyPhiRMS;
}
///_______________________________________________________________________
 Float_t Cluster::GetEnergyRRMS(){
  return m_eRRMS;
}
///_______________________________________________________________________
 Float_t Cluster::GetClusterStartDepth(){
  // Returns the depth where the cluster energy goes above min-I
  Int_t nentries = HitList()->GetEntries();
  Float_t minI = 6.3,depth = 999.;
  Int_t minsys = 2;
  if(!m_clusterStartDepth){
    for(Int_t i=0; i<nentries; i++){
      Float_t E = ((EDeposit*)MCEList()->At(i))->GetEDeposit();
      towerID* tower = (towerID*)HitList()->At(i);
      Int_t sys = tower->system();
      Int_t layer = tower->layer();
      if((E > minI) && (sys <= minsys) && ((sys < minsys)||(layer < depth))){
	depth = layer;
      }
    }
    return(depth);
  }
  else return m_clusterStartDepth;
}
///_______________________________________________________________________
 void Cluster::SetUpCluster(Float_t energy,Float_t ePhi,Float_t eTheta,
		       Float_t R,Float_t phiRMS,Float_t thetaRMS,Float_t RRms){
  // Set the clusters attributes
  m_energy = energy;
  m_energyPhi = ePhi;
  m_energyTheta = eTheta;
  m_eR = R;
  m_energyPhiRMS = phiRMS;
  m_energyThetaRMS = thetaRMS;
  m_eRRMS = RRms;
}
//________________________________________________________________________
 void Cluster::AddHit(towerID* tower,IntObj* index,EDeposit* edep){
  // Add a hit to the cluster
  HitList()->Add(tower);
  IndexList()->Add(index);
  MCEList()->Add(edep);
}
//________________________________________________________________________
 void Cluster::AddParticle(IntObj* part){
  // Add a particle to the cluster
  McList()->Add(part);
}
//________________________________________________________________________
 void Cluster::Spew(FILE* ofile){
  // This is the output function

  fprintf(ofile,"%f %f %f %f %f %f %fn",GetEnergy(),GetEnergyPhi(),
	  GetEnergyTheta(),GetEnergyR(),GetEnergyPhiRMS(),GetEnergyThetaRMS(),
	  GetEnergyRRMS());
	  
  Int_t nentries = HitList()->GetEntries();
  fprintf(ofile,"%i hit(s)n",nentries);

  for(Int_t i=0; i<nentries; i++){
    fprintf(ofile,"%i %x %fn",((IntObj*)m_IndexList[i])->Getnum(),
	    ((towerID*)m_hitList[i])->tag(),
	    ((EDeposit*)m_mcEList[i])->GetEDeposit());
  }

  nentries = McList()->GetEntries();
  fprintf(ofile,"%i particle(s)n",nentries);

  for(Int_t j=0; j<nentries; j++){
    IntObj* part = (IntObj*) McList()->At(j);
    fprintf(ofile,"%in", part->Getnum());
  }
  return;
}