#include "Cluster.h"
#include "McPart.h"
#include "EDeposit.h"
#include "TMath.h"
#include "TObjArray.h"
#include "CalHit.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(TObjArray* ch):m_energy(0),m_energyTheta(0),m_energyPhi(0),m_energyThetaRMS(0),m_energyPhiRMS(0),m_eR(0),m_eRRMS(0),m_clusterStartDepth(0){
  //  constructor when using hits

  m_CalHit = ch;
}
//_________________________________________________________________________
 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();
}
///_______________________________________________________________________
  
 void Cluster::CalcEnergy(){
  // Returns the total energy in the cluster
Float_t Etot = 0.0;
   Int_t nentries = IndexList()->GetEntries();
   for(Int_t i=0; i < nentries; i++){
      Int_t index = ((IntObj*)IndexList()->At(i))->Getnum();
     Etot += ((CalHit*)m_CalHit->At(index))->GetEtot();
   }
   m_energy = Etot;

}
///_______________________________________________________________________
 void Cluster::CalcEnergyTheta(){
  // Returns the energy weighted theta average of the cluster
  Int_t nentries = IndexList()->GetEntries();
  Float_t sumEtheta=0., sumEtheta2=0., sumWidth = 0., sumE=0.;
    for(Int_t i=0; i<nentries; i++){
      Int_t index = ((IntObj*)IndexList()->At(i))->Getnum();
      CalHit* ch = (CalHit*)m_CalHit->At(index);

      Float_t E = ch->GetEtot();
      Float_t ctheta = ch->GetcosTheta();
      sumEtheta += E*ctheta;
      sumEtheta2 += E*ctheta*ctheta;
      Float_t width = ch->GetcosThetaSeg();
      sumWidth += E*width*width;
      sumE += E;
    }
    m_energyTheta = sumEtheta/sumE;
    m_energyThetaRMS = sqrt(sumEtheta2/sumE - 
			    (sumEtheta*sumEtheta)/(sumE*sumE) +
			    sumWidth/sumE/12.);
}
///_______________________________________________________________________
 void Cluster::CalcEnergyPhi(){
  // Returns the energy weighted phi average of the cluster
  Int_t nentries = IndexList()->GetEntries();
  Float_t sumEphi=0., sumEphi2=0., sumWidth=0., sumE = 0., phiRef=0.;
    for(Int_t i=0; i<nentries; i++){

      Int_t index = ((IntObj*)IndexList()->At(i))->Getnum();
      CalHit* ch = (CalHit*)m_CalHit->At(index);
      Float_t E = ch->GetEtot();
      Float_t phi = ch->Getphi();

      if (i==0) {phiRef = phi;}
      Float_t phiDiff = phi - phiRef;

      if (phiDiff > TMath::Pi()) { phiDiff -= 2.*TMath::Pi();}
      else if (phiDiff < -TMath::Pi()) { phiDiff += 2.*TMath::Pi();}

      sumEphi += phiDiff*E;
      sumEphi2 += phiDiff*phiDiff*E;
      Float_t width = ch->GetphiSeg();
      sumWidth += width*width*E;
      sumE += E;
    }
    m_energyPhi = sumEphi/sumE + phiRef;
    m_energyPhiRMS = sqrt(sumEphi2/sumE - (sumEphi*sumEphi)/(sumE*sumE) +
      sumWidth/sumE/12.);
}
///_______________________________________________________________________
 void Cluster::CalcEnergyR(){
  m_eR = 0.;
  m_eRRMS = 0.;
}
 void Cluster::CalcClusterStartDepth(){
  // 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;
     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;
      }
    }
   m_clusterStartDepth = depth;
}
///_______________________________________________________________________
 void Cluster::CalcProperties(){
    CalcEnergy();
    CalcEnergyTheta();
    CalcEnergyPhi();
    CalcClusterStartDepth();
    CalcEnergyR();
};

///_______________________________________________________________________
 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,"%e %e %e %e %e %e %en",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 %en",((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;
}