// ----------------------------------------------------------------------------
// $Id: LCDClusterSystem.cxx,v 1.4 2001/06/21 04:15:50 toshi Exp $
// ----------------------------------------------------------------------------
//
// $Log: LCDClusterSystem.cxx,v $
// Revision 1.4  2001/06/21 04:15:50  toshi
// Changes according to updated LCDCalHit and LCDCalHitUtil.
//
// Revision 1.3  2001/06/20 22:55:01  toshi
// ?????
//
// Revision 1.2  2001/06/20 03:36:47  toshi
// Changes to prevent compilar error on windows.
//
// Revision 1.1  2001/05/10 18:29:03  toshi
// C++ file name convention has been changed from *.C to *.cxx .
//
//

/*
  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 "LCDMcPart.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_calhitutil.SetDetectorParameter(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);
}
//_________________________________________________________________________

//_________________________________________________________________________
 void LCDClusterSystem::SetDetectorParameters(LCDGetParameters* gp) {
  m_calhitutil.SetDetectorParameter(gp);
}
//_________________________________________________________________________

//_________________________________________________________________________
 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;
  Int_t         nparts;
  Double_t      maxE;
  Int_t         max_index;
  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 (m_calhitutil.GetEtot(calhit) < m_ethreshold) continue;//Threshold cut
    nparts = calhit->GetMCEEntries();

    maxE = -100.;  max_index = -1;
    for(j=0 ; j < nparts ; j++){
      partE     = calhit->GetMCEEnergyAt(j);
      partindex = calhit->GetMCEPartAt(j);
      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      = calhit->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++){
      partE     = calhit->GetMCEEnergyAt(j);
      partindex = calhit->GetMCEPartAt(j);
      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       = m_calhitutil.GetEtot(ch);
    theta   = m_calhitutil.GetTheta(ch);
    phi     = TVector2::Phi_mpi_pi(m_calhitutil.GetPhi(ch));
    r       = m_calhitutil.GetR(ch);
    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 = m_calhitutil.GetThetaSeg(ch);
    sumTWidth += t_width*e*t_width;
    sumEphi  += phiDiff*e;
    sumEphi2 += phiDiff*e*phiDiff;
    p_width = m_calhitutil.GetPhiSeg(ch);
    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);

}
//_______________________________________________________________________