//---------------------------------------------------------------------
// File and Version Information:
// 	$Id: $
//
// Description:
//	This object constructs a list of clusters from a list of digis.
//	Online version.
//
// Environment:
//	Software developed for the PANDA Detector at GSI.
//


#include "PndEmcMakeClusterOnline.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 "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"

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


#include <iostream>

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


static Int_t evtCounter = 0;
static Int_t digiCounter = 0;
static Int_t nProg = 0;

PndEmcMakeClusterOnline::PndEmcMakeClusterOnline(Int_t verbose, Bool_t storeclusters):
	FairTask("EmcClusteringTask", verbose),
	fDigiArray(NULL), 
	fClusterArray(NULL), 
	fDigiFunctor(NULL), 
	fDigiEnergyTresholdBarrel(0.), fDigiEnergyTresholdFWD(0.), fDigiEnergyTresholdBWD(0.), fDigiEnergyTresholdShashlyk(0.), 
	fClusterActiveTime(0.),
	fGeoPar(new PndEmcGeoPar()), 
	fRecoPar(new PndEmcRecoPar()), 
	fStoreClusters(storeclusters), 
	fStoreClusterBase(kTRUE),
	fClusterPosParam(),
	fNrOfEvents(0),
	fTimeBetweenDigis(5.)
{
	fClusterPosParam.clear();
}

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

PndEmcMakeClusterOnline::~PndEmcMakeClusterOnline()
{ 
}

// -----   Public method Init   -------------------------------
InitStatus PndEmcMakeClusterOnline::Init() {

	// Get RootManager
	FairRootManager* ioman = FairRootManager::Instance();
	if ( ! ioman )
	{
		cout << "-E- PndEmcMakeClusterOnline::Init: " << "RootManager not instantiated!" << endl;
		return kFATAL;
	}

	/*if(!FairRunAna::Instance()->IsTimeStamp()) {
		cout << "-E- PndEmcMakeClusterOnline::Init: " << "This task can be activated only online" << endl;
		return kFATAL;
	}*/

	// Get nr of events	
	fNrOfEvents = (ioman->GetInTree())->GetEntriesFast();

	// Get input array
	if(FairRunAna::Instance()->IsTimeStamp()) fDigiArray = (TClonesArray*) ioman->GetObject("EmcDigiSorted"); // for timebased, use sorted digis
	else fDigiArray = (TClonesArray*) ioman->GetObject("EmcDigi"); // for eventbased, use "normal" digis
	if (!fDigiArray) {
		cout << "-W- PndEmcMakeClusterOnline::Init: " << "No PndEmcDigi array!" << endl;
		return kERROR;
	}

	// Create and register output array
	fClusterArray = new TClonesArray("PndEmcCluster");
	ioman->Register("EmcClusterTemp","Emc", fClusterArray, fStoreClusters);

	// for subsequent methods we need the "event grouping" as given by the TS buffer 
	// --> fDigiArray becomes an output array. 
	//
	// can be removed once PndEmcCluster object has been rewritten
	if(FairRunAna::Instance()->IsTimeStamp()) {	
		// for subsequent methods we need the "event grouping" as given by the TS buffer --> fDigiArray becomes an output array. 
		fDigiArray = new TClonesArray(fDigiArray->GetClass());		//still needed to activate array 
		ioman->Register("EmcDigiClusterBase", "Emc", fDigiArray, fStoreClusterBase);
		cout << "\n[INFO\t] Running TIMEBASED version." << endl;
	}
	else cout << "\n[INFO\t] Running EVENTBASED version." << endl;

	// searching for time gaps in digi stream
	fDigiFunctor = new TimeGap();

	fGeoPar->InitEmcMapper();  
	PndEmcStructure::Instance();

	// get some parameters from the parbase
	fDigiEnergyTresholdBarrel=fRecoPar->GetEnergyThresholdBarrel();
	fDigiEnergyTresholdFWD=fRecoPar->GetEnergyThresholdFWD();
	fDigiEnergyTresholdBWD=fRecoPar->GetEnergyThresholdBWD();
	fDigiEnergyTresholdShashlyk=fRecoPar->GetEnergyThresholdShashlyk();

	//convert from seconds to nanoseconds...in parfile everything is seconds...here we need nanoseconds
	if (fClusterActiveTime == 0.) fClusterActiveTime=fRecoPar->GetClusterActiveTime() * 1.0e9; 

	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 << "fClusterActiveTime: " << fClusterActiveTime << " ns"<<endl;
	cout << "fTimeBetweenDigis: " << fTimeBetweenDigis << " ns"<<endl;

	cout << fNrOfEvents << " events." << endl;

	cout << "-I- PndEmcMakeClusterOnline: Intialization successful" << endl;
	return kSUCCESS;
}


void PndEmcMakeClusterOnline::Exec(Option_t* opt) {

/*	if (evtCounter == nProg*fNrOfEvents/100) {
		if(FairRunAna::Instance()->IsTimeStamp()) cout << "[INFO\t] " << evtCounter << " timebunches processed." << endl;	
		else cout << "[INFO\t] " << evtCounter << " events processed." << endl;
		nProg += 5;	
	}*/

	if (fVerbose>2){
		fTimer.Start();
	}

	// Reset output arrays and get all digis up to a certain time gap specified by fClusterActiveTime
	if ( ! fClusterArray ) Fatal("Exec", "No Cluster Array");
	fClusterArray->Delete();
	if(FairRunAna::Instance()->IsTimeStamp()) {
		fDigiArray->Delete();
		fDigiArray->AbsorbObjects(FairRootManager::Instance()->GetData("EmcDigiSorted", fDigiFunctor, fClusterActiveTime));
	}	


// --------------------- START OF CLUSTER FINDING ALGORITHM ---------------------------

	Int_t dx = 0;
	Int_t dy = 0;
	Int_t clusterNr = 0;
	Int_t nDigisPassed = 0; // keeps track of the number of digis in this timebunch AND that passed the thresholds, so use as upper limit for loops later on

	Int_t deltaD = 1; // neighbour distance to consider

	Double_t dt = 0;
	Double_t deltaT = 0;

	if(FairRunAna::Instance()->IsTimeStamp()) {
		dt = 0;
		deltaT = fTimeBetweenDigis; // threshold for time between digis in a cluster
	}		

	Int_t nDigis = fDigiArray->GetEntriesFast();

	std::vector<Int_t> neighbours; // make arrays dynamic, as we don't know its size in advance
	std::vector<Int_t> DigiPassed;
	std::vector<Int_t> XPadPassed;
	std::vector<Int_t> YPadPassed;
	std::vector<Double_t> TPassed;
	std::vector<Int_t> tmpArr;
	std::vector<Int_t> isAdded;
	std::vector<Int_t> similarities;

	for (Int_t a=0; a<nDigis; a++) isAdded.push_back(-1); // initialise all entries of isAdded to -1 (needed by algorithm)

	if (fVerbose>2){
		cout<<"DigiList length: "<<nDigis<<endl;
	}

	//loop over all digis to add them to clusters
	
	for (Int_t iDigi=0; iDigi<nDigis; ++iDigi) {

		/* Get a new digi from the digi array */
		PndEmcDigi* theDigi = (PndEmcDigi*) fDigiArray->At(iDigi);

		// thresholds: if energy is below the corresponding threshold, digi is not considered in clustering
		Int_t module=theDigi->GetModule();
		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;

		// fill arrays with digi information for neighbour relationship building
		DigiPassed.push_back(iDigi); // keep track which digis in fDigiArray passed the thresholds, needed to correct for gaps in isAdded[] later on caused by hits being skipped when they are below threshold
		if(FairRunAna::Instance()->IsTimeStamp()) TPassed.push_back(theDigi->GetTimeStamp());
		if (theDigi->GetXPad()<0 && theDigi->GetModule() == 3) XPadPassed.push_back(theDigi->GetXPad()+1); // correct for gaps in FwEndcap
		else XPadPassed.push_back(theDigi->GetXPad());
		if (theDigi->GetYPad()<0 && theDigi->GetModule() == 3) YPadPassed.push_back(theDigi->GetYPad()+1);
		else YPadPassed.push_back(theDigi->GetYPad());
		nDigisPassed++;  
	}

	if (fVerbose>4){
		cout<<"Looped over all digis ("<< nDigisPassed<< " of " << nDigis <<" digis passed)"<<endl;
	}

	if (nDigisPassed>nDigis) {
		Fatal("Exec", "Attempt to process more digis than are present in this timebunch");
	}

//----- First, construct matrix containing information which digis are neighbouring -----//

	for (Int_t d=0; d<nDigisPassed-1; d++) { // check all pairs of digis, so all digis up to the second-to-last one
		Int_t nNeighbours = 0;
		tmpArr.clear();
		for (Int_t e=d+1; e<nDigisPassed; e++) { 
			dx = XPadPassed[e] - XPadPassed[d]; // check euclidian distance between all digis which passed threshold
			dy = YPadPassed[e] - YPadPassed[d]; // this construction ensures all pairs of digis are only checked once
			if(FairRunAna::Instance()->IsTimeStamp()) dt = TMath::Abs(TPassed[e] - TPassed[d]);	// also take into account that non-consecutive digi pairs may differ in time more than dtau ns (timebased only)

			if (dy*dy <= deltaD && dx*dx <= deltaD && dt <= deltaT) { 	// if the distance is small enough, the digis are neighbours (dt = deltaT = 0 for eventbased)
				tmpArr.push_back(e); 		// construct array containing neighbouring digis
				nNeighbours++;
			}
		}
		neighbours.push_back(nNeighbours); // keep track of nr of neighbours
		for (Int_t i=0; i<nNeighbours; i++) {
			neighbours.push_back(tmpArr[i]); 		// so, the stucture of neighbours[] is now [#neighbours_1 hit1 ... hitN #neighbours_2 hit2 ... etc]
		}
	}

	// Primary Clustering
	Int_t k = 0;
	Int_t l = 0;
	Int_t nClusters = 0; // nr of preclusters
	Int_t simLength = 0;
	while (l < neighbours.size() ) {
		if (isAdded[k] < 0) { 		// if it hasn't been added yet
			isAdded[k] = nClusters; // put internal cluster nr in isAdded array
			for (Int_t j=1; j<neighbours[l]+1; j++) {			// for the first neighbouring hit upto #neighbours for this digi
				Int_t m = neighbours[l+j];					// m = hit index
				if (isAdded[m] < 0) isAdded[m] = nClusters; // if hit hasn't been added, put internal cluster nr here
				else if (isAdded[m] != nClusters) {
					similarities.push_back(nClusters);  // if it has, put current cluster nr in this array
					similarities.push_back(isAdded[m]); // together with its cluster nr (could be different, which is why is being stored here for later comparison)
					simLength += 2; // increment simLength by 2, as we just wrote 2 elements to similarities
					}
				}
			nClusters++;
			} 
		else { // if isAdded[k] isn't -1, it means the digi has been added already and isAdded[k] contains its cluster nr
			for (Int_t j=1; j<neighbours[l]+1; j++){
				Int_t m = neighbours[l+j];
				if (isAdded[m]<0) isAdded[m] = isAdded[k]; 	// same as before, only use cluster nr found in isAdded[k]
				else if (isAdded[m] != isAdded[k]) {
					similarities.push_back(isAdded[k]); // similarities[] keeps tracks of which preclusters have to be merged
					similarities.push_back(isAdded[m]); 
					simLength += 2;
					}
				}
			}
		k++;
		l += neighbours[l]+1;
		}

	if (nDigisPassed!=0)
		if (isAdded[nDigisPassed-1]<0) isAdded[nDigisPassed-1] = nClusters++;

    // Secondary clustering
    for (Int_t i=0; i<simLength; i+=2) { // use info from similarities[] to merge preclusters
		for (Int_t m=0; m<nDigisPassed; m++)
			if (isAdded[m] == similarities[i])
				isAdded[m] = similarities[i+1];
    }
    for (Int_t i=0; i<nClusters; i++) {
		Int_t n = 0;
		for (Int_t j=0; j<nDigisPassed; j++)
			if (isAdded[j]==i) n++;
		if (n > 0) {
			for (Int_t j=0; j<nDigisPassed; j++)
				if (isAdded[j]==i) isAdded[j] = clusterNr;
			clusterNr++;
		}
	}

	//--- Finally, use isAdded to merge digis into clusters ---//

	for (Int_t i=0; i<nDigisPassed; i++) {
		if (isAdded[i] < fClusterArray->GetEntriesFast() ) { // if isAdded[i] is smaller than the current #clusters, the cluster to which this digi belongs already exists, and we should add it to this cluster
			PndEmcCluster* cluster= (PndEmcCluster*) fClusterArray->At(isAdded[i]); // set pointer to the cluster we need, as indicated by isAdded[]
			cluster->addDigi(fDigiArray,DigiPassed[i]); // add digi with index as indicated by DigiPassed, so we add the correct one to the cluster
			cluster->AddLink(FairLink("EmcDigi", DigiPassed[i]));
		}
		else { // if not, make a new cluster for this digi
			PndEmcCluster* newCluster = new((*fClusterArray)[fClusterArray->GetEntriesFast()]) PndEmcCluster();
			newCluster->addDigi(fDigiArray, DigiPassed[i]);
			newCluster->SetLink(FairLink("EmcDigi", DigiPassed[i]));
		}
	}

// --------------------- stop---------------------------------------------

	if (evtCounter == nProg*fNrOfEvents/100) {
		if(FairRunAna::Instance()->IsTimeStamp()) cout << "[INFO\t] " << evtCounter << " timebunches processed." << endl;	
		else cout << "[INFO\t] " << evtCounter << " events processed." << endl;
		nProg += 5;	
	}

	evtCounter++;
	
	FinishClusters();

}


void PndEmcMakeClusterOnline::FinishClusters() {

	Int_t nCluster = fClusterArray->GetEntriesFast();
	for (Int_t i=0; i<nCluster; i++) {
		FinishCluster((PndEmcCluster*) (fClusterArray->At(i)));
	}


	if (fVerbose>2){
		fTimer.Stop();
		Double_t rtime = fTimer.RealTime();
		Double_t ctime = fTimer.CpuTime();
		fTimer.Reset();
		cout << "PndEmcMakeClusterOnline, Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
	}	
}

void PndEmcMakeClusterOnline::FinishCluster(PndEmcCluster* cluster) {
	
	std::vector<Int_t> list = cluster->DigiList();

	Double_t total_energy = 0;
	Double_t max_energy = 0;
	Int_t max_energy_idx = 0;

	for(Int_t iDigi=0; iDigi<list.size(); ++iDigi) {
		Int_t idx = list[iDigi];
		PndEmcDigi* thedigi = (PndEmcDigi*) fDigiArray->UncheckedAt(idx);

		total_energy += thedigi->GetEnergy();
		if(thedigi->GetEnergy() > max_energy) {
			max_energy = thedigi->GetEnergy();
			max_energy_idx = idx;
		}
	}

	cluster->SetEnergy(total_energy);

	PndEmcClusterProperties clustProperties(*cluster, fDigiArray);

	TVector3 tmpbumppos = clustProperties.Where(fRecoPar->GetEmcClusterPosMethod(), fClusterPosParam);
	cluster->SetPosition(tmpbumppos);
	
	Double_t distanceToIP = tmpbumppos.Mag()/100; // distance to IP in m (default position is in cm)
	Double_t flightTime = distanceToIP/.299792458; // flight time in ns = distance to IP / c

	cluster->SetTimeStamp(static_cast<PndEmcDigi*>(fDigiArray->UncheckedAt(max_energy_idx))->GetTimeStamp()-flightTime); // correct cluster time for photon flight time

	PndEmcXClMoments xClMoments(*cluster, fDigiArray);
	cluster->SetZ20(xClMoments.AbsZernikeMoment(2, 0, 15));
	cluster->SetZ53(xClMoments.AbsZernikeMoment(5, 3, 15));
	cluster->SetLatMom(xClMoments.Lat());

	//-------------------stop-------------------------------------------------------------
}


void PndEmcMakeClusterOnline::FinishTask() {
	if(fVerbose>2) {
		std::cout << "Processed in total " << digiCounter << " digis within " << evtCounter << " events." << std::endl;
	}
}

void PndEmcMakeClusterOnline::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 reconstruction parameter container
	fRecoPar = (PndEmcRecoPar*) db->getContainer("PndEmcRecoPar");
}


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



ClassImp(PndEmcMakeClusterOnline)