//---------------------------------------------------------------------
// File and Version Information:
// 	$Id: $
//
// Description:
//	This object constructs a list of clusters from a list of digis.
//      This module use a simple algorithm.
//      Digis which are in a Cluster are add to this cluster.
//      If a Digi is in more than one cluster, the clusters are merged.
//      Digis are added to the cluster if they have an energy above
//      _digiEnergyTreshold->value().
//
// Environment:
//	Software developed for the PANDA Detector at GSI.
//
// Author List:
//	Jan Zhong            //------------------------------------------------------------------------

#include "PndEmcMakeCluster.h"

#include "PndEmcClusterProperties.h"
#include "PndEmcXClMoments.h"
#include "PndEmcStructure.h"
#include "PndEmcMapper.h"
#include "PndEmcGeoPar.h"
#include "PndEmcDigiPar.h"
#include "PndEmcRecoPar.h"
#include "PndEmcCluster.h"
#include "PndEmcDigi.h"
#include "PndMCTrack.h"

#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"

#include "TClonesArray.h"
#include "TStopwatch.h"
#include "TROOT.h"

#include <iostream>
//#include <vector>
//#include <set>

using std::cout;
using std::endl;

Int_t PndEmcMakeCluster::fEventCounter=0;

PndEmcMakeCluster::PndEmcMakeCluster(Int_t verbose, Bool_t storeclusters) 
{
	fVerbose=verbose;
	fStoreClusters=storeclusters;
}

//--------------
// Destructor --
//--------------

PndEmcMakeCluster::~PndEmcMakeCluster()
{ 
}

// -----   Public method Init   -------------------------------
InitStatus PndEmcMakeCluster::Init() {
	
  	// Get RootManager
	FairRootManager* ioman = FairRootManager::Instance();
	if ( ! ioman )
	{
		cout << "-E- PndEmcMakeCluster::Init: "
		<< "RootManager not instantiated!" << endl;
		return kFATAL;
	}
	
	// Get input array
	fDigiArray = (TClonesArray*) ioman->GetObject("EmcDigi");
	if ( ! fDigiArray ) {
		cout << "-W- PndEmcMakeCluster::Init: "
		<< "No PndEmcDigi array!" << endl;
		return kERROR;
	}
	
	// Get input array
	fHitArray = (TClonesArray*) ioman->GetObject("EmcHit");
	if ( ! fHitArray ) {
		cout << "-W- PndEmcMakeCluster::Init: "
		<< "No PndEmcHit array! Needed for MC Truth" << endl;
//		return kERROR;
	}

	// Get input array
	fMCTrackArray = (TClonesArray*) ioman->GetObject("MCTrack");
	if ( ! fMCTrackArray ) {
		cout << "-W- PndEmcMakeCluster::Init: "
		<< "No MCTrack array! Needed for MC Truth" << endl;
//		return kERROR;
	}

	// Create and register output array
	fClusterArray = new TClonesArray("PndEmcCluster");
	
	ioman->Register("EmcCluster","Emc",fClusterArray,fStoreClusters);
	
	fGeoPar->InitEmcMapper();  
	PndEmcStructure::Instance();

	fDigiEnergyTresholdBarrel=fRecoPar->GetEnergyThresholdBarrel();
	fDigiEnergyTresholdFWD=fRecoPar->GetEnergyThresholdFWD();
	fDigiEnergyTresholdBWD=fRecoPar->GetEnergyThresholdBWD();
	fDigiEnergyTresholdShashlyk=fRecoPar->GetEnergyThresholdShashlyk();
	
	cout<<"PndEmcMakeCluster::fDigiEnergyTresholdBarrel: "<<fDigiEnergyTresholdBarrel<<endl;
	cout<<"PndEmcMakeCluster::fDigiEnergyTresholdFWD: "<<fDigiEnergyTresholdFWD<<endl;
	cout<<"PndEmcMakeCluster::fDigiEnergyTresholdBWD: "<<fDigiEnergyTresholdBWD<<endl;
	cout<<"PndEmcMakeCluster::fDigiEnergyTresholdShashlyk: "<<fDigiEnergyTresholdShashlyk<<endl;

	if (!strcmp(fRecoPar->GetEmcClusterPosMethod(),"lilo"))
	{
		cout<<"Lilo cluster position method"<<endl;
		fClusterPosParam.push_back(fRecoPar->GetOffsetParmA());
		fClusterPosParam.push_back(fRecoPar->GetOffsetParmB());
		fClusterPosParam.push_back(fRecoPar->GetOffsetParmC());
	}	

	cout << "-I- PndEmcMakeCluster: Intialization successfull" << endl;

	return kSUCCESS;
}


void PndEmcMakeCluster::Exec(Option_t* opt) 
{
	TStopwatch timer;
	if (fVerbose>2){
		timer.Start();
	}
	// Reset output array
	if ( ! fClusterArray ) Fatal("Exec", "No Cluster Array");
	fClusterArray->Delete();

	Int_t nDigis = fDigiArray->GetEntriesFast();
	if (fVerbose>2){
		cout<<"DigiList length "<<nDigis<<endl;
	}
	
	Double_t totalDigiEnergy=0;
	//loop to build Cluster
	for (Int_t iDigi=0; iDigi<nDigis; iDigi++)
	{
		PndEmcDigi* theDigi = (PndEmcDigi*) fDigiArray->At(iDigi);
		
		Int_t module=theDigi->GetModule();
		
		// In the following lines there is separate threshold for forward endcup
		// and the same for barrel and backward and shashlyk EMC, but for last 2 threshold probably should be different 
		if ((module==1||module==2)&&(theDigi->GetEnergy()< fDigiEnergyTresholdBarrel)) continue;
		if ((module==3)&&(theDigi->GetEnergy()< fDigiEnergyTresholdFWD)) continue;
		if ((module==4)&&(theDigi->GetEnergy()< fDigiEnergyTresholdBWD)) continue;
		if ((module==5)&&(theDigi->GetEnergy()< fDigiEnergyTresholdShashlyk)) continue;
		
			bool isAdded = false;
		
			Int_t clustmarker=0;
			
			Int_t clustLength=fClusterArray->GetEntriesFast();
			
			for(Int_t i=0;i<clustLength;i++)
			{
				PndEmcCluster* cluster=(PndEmcCluster*) fClusterArray->At(i);
				if(cluster->isInCluster(theDigi, fDigiArray))
				{
					if(!isAdded)
					{
						clustmarker=i;
						isAdded=true;
						cluster->addDigi(fDigiArray, iDigi);
						cluster->AddLink(FairLink("EmcDigi", iDigi));
					}
					else
					{
						PndEmcCluster* clust_clustmarker=(PndEmcCluster*) fClusterArray->At(clustmarker);
						PndEmcCluster* clust_i=(PndEmcCluster*) fClusterArray->At(i);
						clust_clustmarker->addCluster(clust_i, fDigiArray);
						fClusterArray->RemoveAt(i);
						fClusterArray->Compress();
						clustLength--;
						i--;
					}
				}
			}
	
			if (!isAdded)
			{
				PndEmcCluster* newcluster = new((*fClusterArray)[clustLength]) PndEmcCluster();
				newcluster->addDigi(fDigiArray, iDigi);
				newcluster->SetLink(FairLink("EmcDigi", iDigi));
			}
			
			totalDigiEnergy+=theDigi->GetEnergy();
			
	}
	
	// At that moment internal state fEnergy and fWhere of Clusters are not initialized, the following does initialisation
	Int_t nCluster = fClusterArray->GetEntriesFast();
	for (Int_t i=0; i<nCluster; i++)
	{
		PndEmcCluster *tmpclust = (PndEmcCluster*) fClusterArray->At(i);
		PndEmcClusterProperties clustProperties(*tmpclust, fDigiArray);

		tmpclust->SetEnergy(clustProperties.Energy());
		TVector3 tmpbumppos = clustProperties.Where(fRecoPar->GetEmcClusterPosMethod(), fClusterPosParam);
		tmpclust->SetPosition(tmpbumppos);
		PndEmcXClMoments xClMoments(*tmpclust, fDigiArray);
		tmpclust->SetZ20(xClMoments.AbsZernikeMoment(2, 0, 15));
		tmpclust->SetZ53(xClMoments.AbsZernikeMoment(5, 3, 15));
		tmpclust->SetLatMom(xClMoments.Lat());
		tmpclust->fMcList.clear();
		if(fHitArray && fMCTrackArray){
			// BS: this is a first order approximation only !!!!

			std::vector <Int_t> newlist;
			newlist.clear();
			for(Int_t j=0; j<tmpclust->fDigiList.size(); j++){
				PndEmcDigi*m;
				m=(PndEmcDigi*)fDigiArray->At(tmpclust->fDigiList[j]);

				Int_t inx;
				inx=m->GetHitIndex();

				if(inx>=0){
					const std::vector <Int_t> &tmplist=((PndEmcHit*) fHitArray->At(inx))->GetMcList();

					// I copy the complete list instead of only the highest energy particle
					// highest energy might be a problem for Hits which belong to two real clusters?
					for(Int_t k=0; k<tmplist.size(); k++){
						newlist.push_back(tmplist[k]);
					}
				}
			}
			// if( fMCTrackArray) checked above already
	        cleansortmclist(newlist,fMCTrackArray);
			tmpclust->fMcList=newlist;
		}
	}
	
	fEventCounter++;
	if (fVerbose>0)
		cout<<"PndEmcMakeCluster, event: "<<fEventCounter<<endl;
	
	if (fVerbose>2){
		timer.Stop();
		Double_t rtime = timer.RealTime();
		Double_t ctime = timer.CpuTime();
		cout << "PndEmcMakeCluster, Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
	}	
}

// Helper function, does not depend on class, identical to the one in PndEmcHitProducer
void PndEmcMakeCluster::cleansortmclist( std::vector <Int_t> &newlist,TClonesArray* mcTrackArray)
{
	std::vector <Int_t> tmplist;
	// Sort list...
	std::sort( newlist.begin(), newlist.end());
	// and copy every id only once (even so it might be in the list several times)
	std::unique_copy( newlist.begin(), newlist.end(), std::back_inserter( tmplist ) );

	// Now check if mother or (grand)^x-mother are already in the list
	// (which means i am a secondary)... if so, remove myself
	for(Int_t j=tmplist.size()-1; j>=0; j--){
		bool flag;
		PndMCTrack *pt;
		pt=((PndMCTrack*)mcTrackArray->At(tmplist[j]));
		if(pt->GetMotherID()<0) continue;
		flag=false;
		while(!flag){
			Int_t id;
			id=pt->GetMotherID();
			if(id<0) break;
			pt=(PndMCTrack*)mcTrackArray->At(id);

			for(Int_t k=j-1; k>=0; k--){
				if(tmplist[k]==id){
					tmplist.erase(tmplist.begin()+j);
					flag=true;
					break;
				}
			}
		}
	}
	newlist=tmplist;
}

void PndEmcMakeCluster::SetParContainers() {

  // Get run and runtime database
  FairRun* run = FairRun::Instance();
  if ( ! run ) Fatal("SetParContainers", "No analysis run");

  FairRuntimeDb* db = run->GetRuntimeDb();
  if ( ! db ) Fatal("SetParContainers", "No runtime database");

  // Get Emc geometry parameter container
  fGeoPar = (PndEmcGeoPar*) db->getContainer("PndEmcGeoPar");
  
  // Get Emc digitisation parameter container
  fDigiPar = (PndEmcDigiPar*) db->getContainer("PndEmcDigiPar");
 
  // Get Emc reconstruction parameter container
  fRecoPar = (PndEmcRecoPar*) db->getContainer("PndEmcRecoPar");
}

void PndEmcMakeCluster::SetStorageOfData(Bool_t val)
{
  fStoreClusters=val;
  return;
}
  

ClassImp(PndEmcMakeCluster)