/*
  ClusterSystem

  The ClusterSystem class is used to extract cluster information
  from the GISMO ascii output file and use it to fill the Cluster class.


  V 1.0 XX-XXXXX-XXXX XXXXXXX    Program creation.
  V 2.0 21-April-2000 M.Iwasaki  Make and fill 3 kinds of clusters: 
                                 EM, HAD and MUClusters separately.
*/

#include <iostream.h>

#include "LCDClusterSystem.h"

#include "LCDCalHit.h"
#include "LCDMcPart.h"
#include "LCDEDeposit.h"

ClassImp(LCDClusterSystem)

//_________________________________________________________________________
 LCDClusterSystem::LCDClusterSystem() {
  // Default constructor
  Init();
}
//_________________________________________________________________________

//_________________________________________________________________________
 LCDClusterSystem::LCDClusterSystem(LCDGetParameters* gp) {
  // Default constructor
  Init();
  SetDetectorParameters(gp);
}
//_________________________________________________________________________

//______________________________________________________________________
 LCDClusterSystem::~LCDClusterSystem(){
  // Cleans up loose ends. There aren't any right now.
  delete m_reference_lm;
  delete m_reference_em_b;
  delete m_reference_em_e;
  delete m_reference_hd_b;
  delete m_reference_hd_e;
  delete m_reference_mu_b;
  delete m_reference_mu_e;
}
//_______________________________________________________________________

//_________________________________________________________________________
 void LCDClusterSystem::Init() {
  m_parameters = 0;
  m_ethreshold = 10e-3;//10MeV
  m_reference_em_b  = new TMap(17,0);
  m_reference_em_e  = new TMap(17,0);
  m_reference_hd_b  = new TMap(17,0);
  m_reference_hd_e  = new TMap(17,0);
  m_reference_mu_b  = new TMap(17,0);
  m_reference_mu_e  = new TMap(17,0);
  m_reference_lm    = new TMap(17,0);
  LCDCalHit::energyScale = m_hitEscale;
  LCDCalHit::ThetaSeg = m_hitThetaSeg;
  LCDCalHit::phiSeg = m_hitphiSeg;
  LCDCalHit::barrelThick = m_barrelThick;
  LCDCalHit::endcapThick = m_endcapThick;
  LCDCalHit::barrelInner = m_barrelInner;
  LCDCalHit::endcapInner = m_endcapInner;
}
//_________________________________________________________________________

//_________________________________________________________________________
 void LCDClusterSystem::SetDetectorParameters(LCDGetParameters* gp) {
  m_parameters = gp;
  m_hitEscale[0] = gp->GetEMEnergyScale();
  m_hitEscale[1] = gp->GetHADEnergyScale();
  m_hitEscale[2] = gp->GetMUEnergyScale();
  m_hitEscale[3] = gp->GetLUMEnergyScale();

  m_hitThetaSeg[0] = gp->GetEMThetaSeg();
  m_hitThetaSeg[1] = gp->GetHADThetaSeg();
  m_hitThetaSeg[2] = gp->GetMUThetaSeg();
  m_hitThetaSeg[3] = gp->GetLUMThetaSeg();

  m_hitphiSeg[0] = gp->GetEMphiSeg();
  m_hitphiSeg[1] = gp->GetHADphiSeg();
  m_hitphiSeg[2] = gp->GetMUphiSeg();
  m_hitphiSeg[3] = gp->GetLUMphiSeg(); 

  m_barrelThick[0] = gp->GetEMBarrelThick();
  m_barrelThick[1] = gp->GetHADBarrelThick();
  m_barrelThick[2] = gp->GetMUBarrelThick();

  m_endcapThick[0] = gp->GetEMEndcapThick();
  m_endcapThick[1] = gp->GetHADEndcapThick();
  m_endcapThick[2] = gp->GetMUEndcapThick();
  m_endcapThick[3] = gp->GetLUMEndcapThick();

  m_barrelInner[0] = gp->GetEMBarrelInnerR();
  m_barrelInner[1] = gp->GetHADBarrelInnerR();
  m_barrelInner[2] = gp->GetMUBarrelInnerR(); 

  m_endcapInner[0] = gp->GetEMEndcapInnerZ();
  m_endcapInner[1] = gp->GetHADEndcapInnerZ();
  m_endcapInner[2] = gp->GetMUEndcapInnerZ();
  m_endcapInner[3] = gp->GetLUMEndcapInnerZ();

}
//_________________________________________________________________________

//_________________________________________________________________________
 void LCDClusterSystem::FillClustersCheat(TObjArray* cal_digi, 
					 TClonesArray* Cluster_list,
					 TObjArray* mclist){

  // Clustering using Cheat method...

  Int_t nhits = cal_digi->GetEntries();

  LCDCalHit*    calhit;
  Int_t         i,j;
  LCDEDeposit*  list;
  Int_t         nparts;
  Double_t      maxE;
  Int_t         max_index;
  LCDEDeposit*  edep;
  Double_t      partE;
  Int_t         partindex;
  LCDMcPart*    part;
  Double_t      dep;
  Int_t         systemid;
  LCDCluster*   kal;
  TMap*         ref;

  Int_t        nkal=0;
  for(i=0 ; i < nhits ; i++){
    calhit = (LCDCalHit*)cal_digi->At(i);
    if (calhit->GetEtot() < m_ethreshold) continue;//Threshold cut
    list   = calhit->GetMCE();
    nparts = calhit->GetMCEEntries();

    maxE = -100.;  max_index = -1;
    for(j=0 ; j < nparts ; j++){
      edep      = &list[j];
      partE     = edep->GetEDeposit();
      partindex = edep->GetPart();
      if (partE > maxE) {
	maxE = partE;  max_index = partindex;
      }
    }

    if (max_index == -1) continue;

    part = (LCDMcPart*)mclist->At(max_index);

    systemid = 
      ((LCDtowerID*)((LCDCalHit*)cal_digi->At(i))->GetTower())->GetSystem();

    dep      = ((LCDCalHit*)cal_digi->At(i))->GetEtotRaw();

    switch(systemid) {
    case  5://LUM now kill LUM
      ref = m_reference_lm;
      break;
    case 10://EM Barrel
      ref = m_reference_em_b;
      break;
    case 11://EM Endcap
      ref = m_reference_em_e;
      break;
    case 12://HAD Barrel
      ref = m_reference_hd_b;
      break;
    case 13://HAD Endcap
      ref = m_reference_hd_e;
      break;
    case 14://MU Barrel
      ref = m_reference_mu_b;
      break;
    case 15://MU Endcap
      ref = m_reference_mu_e;
      break;
    default:
      continue;
      break;
    }

    if(ref->FindObject(part)){// Does cluster exist?

      kal = (LCDCluster*)(ref->GetValue(part));
      kal->AddCalHit(i);

    } else {

      kal = new((*Cluster_list)[nkal++]) LCDCluster();
      kal->AddCalHit(i);
      ref->Add(part,kal);

      switch(systemid) {
      case 5: //LUM
	kal->SetSystem(3);
	kal->SetBarEnd(1);
	break;
      case 10: //EM?
	kal->SetSystem(0);
	kal->SetBarEnd(0);
	break;
      case 11: //EM?
	kal->SetSystem(0);
	kal->SetBarEnd(1);
	break;
      case 12: //HAD?
	kal->SetSystem(1);
	kal->SetBarEnd(0);
	break;
      case 13://HAD?
	kal->SetSystem(1);
	kal->SetBarEnd(1);
	break;
      case 14://MU?
	kal->SetSystem(2);
	kal->SetBarEnd(0);
	break;
      case 15://MU?
	kal->SetSystem(2);
	kal->SetBarEnd(1);
	break;
      }

    }

    for(j=0 ; j < nparts ; j++){
      edep      = &list[j];
      partE     = edep->GetEDeposit();
      partindex = edep->GetPart();
      kal->AddParticle(partindex,partE);
    }

  }

  //Calc...
  for (Int_t ikal=0 ; ikal < nkal ; ikal++) {
    kal=(LCDCluster*)Cluster_list->At(ikal);
    CalcProperties(cal_digi,kal);
  }

  m_reference_lm->Clear();
  m_reference_em_b->Clear();
  m_reference_em_e->Clear();
  m_reference_hd_b->Clear();
  m_reference_hd_e->Clear();
  m_reference_mu_b->Clear();
  m_reference_mu_e->Clear();

}
//_______________________________________________________________________

//_______________________________________________________________________
 void LCDClusterSystem::Doit(LCDEvent* event){
  // This function executes cluster cheat system

  FillClustersCheat(event->CAL(),event->ClusterLst(),event->MCparticles());
}
//_______________________________________________________________________

//_______________________________________________________________________
 void LCDClusterSystem::Cleanup(LCDEvent* event){
  // Cleans up loose ends. There aren't any right now.
  if (event)
    event->ClusterLst()->Delete();
}
//_______________________________________________________________________

//_______________________________________________________________________
 void LCDClusterSystem::CalcProperties(TObjArray* cal_digi,LCDCluster* kal){
  // Returns the total energy in the cluster
  Double_t sumE=0.;
  Double_t sumEtheta=0.;
  Double_t sumEtheta2=0.;
  Double_t sumTWidth = 0.;
  Double_t sumEphi=0.;
  Double_t sumEphi2=0.;
  Double_t sumPWidth=0.;
  Double_t sumEr  = 0.;
  Double_t sumEr2 = 0.;
  Double_t e;
  Double_t theta;
  Double_t phi;
  Double_t r;
  Double_t t_width;
  Double_t p_width;
  Double_t phiRef=0.;
  Double_t phiDiff;
  Double_t rRef  = 0.;
  Double_t rDiff;
  //Double_t minI = 6.3;
  Double_t depth = 999.;
  //Int_t    minsys = 2;
  Int_t    sys;
  Int_t    layer;

  Int_t    nhit   = kal->GetCalHitEntries();
  LCDCalHit* ch;
  for(Int_t i=0; i < nhit ; i++){
    ch = (LCDCalHit*) cal_digi->At(kal->GetCalHitAt(i));
    
    e       = ch->GetEtot();
    theta   = ch->GetTheta();
    phi     = TVector2::Phi_mpi_pi(ch->GetPhi());
    r       = ch->GetR();
    if (i==0) {phiRef = phi;}
    if (i==0) {rRef   = r;  }
    phiDiff = TVector2::Phi_mpi_pi(phi - phiRef);
    rDiff   = r - rRef;

    sumE += e;
    sumEtheta  += theta*e;
    sumEtheta2 += theta*e*theta;
    t_width = ch->GetThetaSeg();
    sumTWidth += t_width*e*t_width;
    sumEphi  += phiDiff*e;
    sumEphi2 += phiDiff*e*phiDiff;
    p_width = ch->GetPhiSeg();
    sumPWidth += p_width*e*p_width;
    sumEr +=rDiff*e;
    sumEr2+=rDiff*e*rDiff;

    LCDtowerID* tower = ch->GetTower();
    sys = tower->GetSystem();
    layer = tower->GetLayer();
    if(e > m_ethreshold && layer < depth ){
      depth = layer;
    }

  }
  kal->SetEnergy(sumE);
  kal->SetEnergyTheta(sumEtheta/sumE);
  kal->SetEnergyThetaRMS(TMath::Sqrt(sumEtheta2/sumE
                                     - (sumEtheta*sumEtheta)/(sumE*sumE)
                                     + sumTWidth/sumE/12.));
  kal->SetEnergyPhi(TVector2::Phi_mpi_pi(sumEphi/sumE + phiRef));
  kal->SetEnergyPhiRMS(TMath::Sqrt(sumEphi2/sumE 
                                   - (sumEphi*sumEphi)/(sumE*sumE)
                                   + sumPWidth/sumE/12.));
  if (sumE != 0.) {
    kal->SetEnergyR(sumEr/sumE + rRef);
    kal->SetEnergyRRMS(TMath::Sqrt(sumEr2/sumE-(sumEr/sumE)*(sumEr/sumE)));
  }
  if (sumE == 0.) {
    kal->SetEnergyR(   10000.0);
    kal->SetEnergyRRMS(10000.0);
  }
  kal->SetClusterStartDepth(depth);

}
//_______________________________________________________________________