//---------------------------------------------------------------------
// File and Version Information:
// 	$Id: $
//
// Description:
//	This object constructs a list of clusters from a list of digis.
// 	+ Implementation of online clustering algorithm
// 	+ Mapping of digis into virtual DC (more closely resembles real readout chain)
//
// Environment:
//	Software developed for the PANDA Detector at GSI.
//
// Note:
//	Because of mapping, probably much slower than other cluster finding programs. Recommended to use for testing of online trigger only.


#include "PndEmcDistributedClustering.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>
#include <vector>

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


static Int_t evtCounter = 0;
static Int_t digiCounter = 0;
static Int_t nProg = 0;
static Double_t DCtotRtime = 0;
static Double_t DCtotCtime = 0;
static Double_t CNtotRtime = 0;
static Double_t CNtotCtime = 0;

PndEmcDistributedClustering::PndEmcDistributedClustering(Int_t verbose, Bool_t storeclusters):
	FairTask("EmcClusteringTask", verbose),
	fDigiArray(NULL), 
	fPreclusterArray(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 --
//--------------

PndEmcDistributedClustering::~PndEmcDistributedClustering()
{ 
}

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

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

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

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

	// Get input array
	fDigiArray = (TClonesArray*) ioman->GetObject("EmcDigiSorted");
	if (!fDigiArray) {
		cout << "-W- PndEmcDistributedClustering::Init: " << "No PndEmcDigi array!" << endl;
		return kERROR;
	}

	// Create and register output array
	fPreclusterArray = new TClonesArray("PndEmcCluster");
	ioman->Register("EmcPrecluster","Emc", fPreclusterArray, fStoreClusters);
	fClusterArray = new TClonesArray("PndEmcCluster");
	ioman->Register("EmcCluster","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
	fDigiArray = new TClonesArray(fDigiArray->GetClass());		//still needed to activate array 
	ioman->Register("EmcDigiClusterBase", "Emc", fDigiArray, fStoreClusterBase);

	// 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- PndEmcDistributedClustering: Intialization successful" << endl;
	return kSUCCESS;
}


void PndEmcDistributedClustering::Exec(Option_t* opt) {

	if (evtCounter == nProg*fNrOfEvents/100) {
		cout << "[INFO\t] " << evtCounter << " timebunches processed." << endl;	
		nProg += 5;	
	}

/*	if (fVerbose>3){
		fTimer.Start();
	}*/

	// Reset output arrays
	if ( ! fClusterArray ) Fatal("Exec", "No Cluster Array");
	fClusterArray->Delete();
	fPreclusterArray->Delete();
	fDigiArray->Delete();

	
	// Get all digis up to a certain time gap specified by fClusterActiveTime
	fDigiArray->AbsorbObjects(FairRootManager::Instance()->GetData("EmcDigiSorted", fDigiFunctor, fClusterActiveTime));

	Int_t nDigis = fDigiArray->GetEntriesFast();

	digiCounter += nDigis;

	// -----   Timer for preclustering  -----------------------------------------
	TStopwatch DCtimer;
	DCtimer.Start();

// --------------------- START SEARCHING FOR PRECLUSTERS ---------------------------

	Int_t dx = 0;
	Int_t dy = 0;
	Int_t nPreclusters = 0; 
	Int_t nDigisPassed = 0; // #digis that passed thresholds
	Int_t nDCdigis = 0;
	Int_t deltaD = 1; // neighbour distance to consider

	Double_t deltaT = fTimeBetweenDigis; // threshold for time between digis in a cluster
//	Double_t deltaT = fClusterActiveTime; // use dtau as threshold for digis inside a cluster as well	
	Double_t dt = 0;

	if (fVerbose>2){
		cout<<"#digis: "<<nDigis<<endl;
		cout << "Per DC: ";
	}

	std::vector<Int_t> neighbours; // make neighbours 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<Int_t> tmpArr;
	std::vector<Int_t> isAdded;
	std::vector<Int_t> similarities;
	std::vector<Double_t> TPassed;
	std::vector<Bool_t> tagArray; 

	for (Int_t iDC=0; iDC<159; iDC++) { // loop over all virtual DCs

	//----- Tag which digis are in the current virtual DC -----//

		tagArray.clear();
		nDCdigis = 0;

		for (Int_t i=0; i<nDigis; i++) {
			PndEmcDigi *myDigi=(PndEmcDigi*)fDigiArray->At(i);
			if (myDigi->GetDCnumber() == iDC) {
				tagArray.push_back(true);
				nDCdigis++;
			}
			else tagArray.push_back(false);
		}

		nDigisPassed = 0; 

		if (nDCdigis == 0) continue; // skip empty DCs

		if (fVerbose>2) cout << nDCdigis << " " << std::flush;

	// reset vectors
		neighbours.clear();
		DigiPassed.clear();
		XPadPassed.clear();
		YPadPassed.clear();
		TPassed.clear();
		isAdded.clear();
		similarities.clear();

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

	//----- Loop over all digis to add them to preclusters -----//
	
		for (Int_t iDigi=0; iDigi<nDigis; ++iDigi) {

			if (tagArray[iDigi] == false) continue; // only process digis that are in the current virtual DC

			/* 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
			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 " << nDCdigis <<" digis passed)"<<endl;
		}

		if (nDigisPassed>nDCdigis) {
			Fatal("Exec", "Attempt to process more digis than are present in this virtual data concentrator");
		}

		if (nDigisPassed==0) continue; // skip DCs that only have digis with energies below threshold

	//----- 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
				dt = TMath::Abs(TPassed[e] - TPassed[d]);	// also take into account that non-consecutive digi pairs may differ in time more than dtau ns

				if (dy*dy <= deltaD && dx*dx <= deltaD && dt <= deltaT) { 	// if the distance is small enough, the digis are neighbours
					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 (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] = nPreclusters; // info in isAdded: digi "index" belongs to precluster "isAdded[index]"
				nPreclusters++;									  // cluster nr, don't reset it while searching for preclusters in a timebunch
			}
		}

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

		for (Int_t i=0; i<nDigisPassed; i++) { 
			if (isAdded[i] < 0) continue;
			if (isAdded[i] < fPreclusterArray->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* precluster= (PndEmcCluster*) fPreclusterArray->At(isAdded[i]); // set pointer to the cluster we need, as indicated by isAdded[]
				precluster->addDigi(fDigiArray,DigiPassed[i]); // add digi with index as indicated by DigiPassed, so we add the correct one to the cluster
				precluster->AddLink(FairLink("EmcDigi", DigiPassed[i]));
			}
			else { // if not, make a new cluster for this digi
				PndEmcCluster* newPrecluster = new((*fPreclusterArray)[fPreclusterArray->GetEntriesFast()]) PndEmcCluster();
				newPrecluster->addDigi(fDigiArray, DigiPassed[i]);
				newPrecluster->SetLink(FairLink("EmcDigi", DigiPassed[i]));
			}
		}

	}

	if (fVerbose>2) cout << endl;

// ----- FINISHED BUILDING PRECLUSTERS -----------------------------------

	if (fVerbose>3){
		cout<<"Finished digi assignment: "<<nPreclusters<<" preclusters identified ("<<std::flush;
	}	
	FinishPreclusters();

	DCtimer.Stop();
	Double_t DCrtime = DCtimer.RealTime();
	Double_t DCctime = DCtimer.CpuTime();
	DCtimer.Reset();

	DCtotRtime += DCrtime; // keep track how much time was spent on preclustering
	DCtotCtime += DCctime;

	// -----   Timer for clustering  -------------------------------------
	TStopwatch CNtimer;
	CNtimer.Start();

// ----- START MERGING PRECLUSTERS INTO 'REAL' CLUSTERS ------------------

	Int_t clusterNr = 0; // #preclusters

	std::vector<Int_t> neighbours2; // make neighbours2 dynamic, as we don't know its size in advance
	std::vector<Int_t> tmpArr2; 
	std::vector<Int_t> isAdded2; 
	std::vector<Int_t> similarities2; 

	Int_t nPres = fPreclusterArray->GetEntriesFast();

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

	if (fVerbose>2){
		cout<<"#preclusters: "<<nPres<<endl;
	}

	//loop over all preclusters to add them to clusters
	
	for (Int_t iClus=0; iClus<nPres-1; ++iClus) {

		Int_t nNeighbours = 0; // reset nNeighbours
		tmpArr2.clear();

		/* Get a new precluster from the precluster array */
		PndEmcCluster* preclus1 = (PndEmcCluster*) fPreclusterArray->At(iClus);

		for (Int_t j=iClus+1; j<nPres; j++) {

			PndEmcCluster* preclus2 = (PndEmcCluster*) fPreclusterArray->At(j);

			if (fVerbose>3){
				cout<<"Distance between preclusters " << iClus << " and " << j << ": "<<std::flush;
				cout<<preclus1->DistanceToCentre(preclus2->where())<<endl;
			}

		// tag which preclusters should be merged
//							spatial distance																							 time difference			(use the same as for digis)
//									|																											|								|
//									V																											V								V
			if ( preclus1->DistanceToCentre(preclus2->where()) <= (preclus1->radius()+preclus2->radius()) && TMath::Abs(preclus1->GetTimeStamp()-preclus2->GetTimeStamp()) < deltaT ) {
				tmpArr2.push_back(j);
				nNeighbours++;
			}
		}

	//--- First, construct vector containing information which preclusters are neighbouring ---//

		neighbours2.push_back(nNeighbours); 		// keep track of nr of neighbours
		for (Int_t i=0; i<nNeighbours; i++) {
			neighbours2.push_back(tmpArr2[i]);	// so, the stucture of neighbours2[] 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 < neighbours2.size() ) {
		if (isAdded2[k] < 0) { 		// if it hasn't been added yet
			isAdded2[k] = nClusters; // put cluster nr in isAdded2 array
			for (Int_t j=1; j<neighbours2[l]+1; j++) {			// for the first neighbouring hit upto #neighbours for this digi
				Int_t m = neighbours2[l+j];					// m = hit index
				if (isAdded2[m] < 0) isAdded2[m] = nClusters; // if hit hasn't been added, put cluster nr here
				else if (isAdded2[m] != nClusters) {
						similarities2.push_back(nClusters);  // if it has, put current cluster nr in this array
						similarities2.push_back(isAdded2[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 isAdded2[k] isn't -1, it means the digi has been added already and isAdded2[k] contains its cluster nr
			for (Int_t j=1; j<neighbours2[l]+1; j++){
				Int_t m = neighbours2[l+j];
				if (isAdded2[m]<0) isAdded2[m] = isAdded2[k]; 	// same as before, only use cluster nr found in isAdded2[k]
				else if (isAdded2[m] != isAdded2[k]) {
						similarities2.push_back(isAdded2[k]); // similarities[] keeps tracks of which preclusters have to be merged
						similarities2.push_back(isAdded2[m]); 
						simLength += 2;
				}
			}
		}
		k++;
		l += neighbours2[l]+1;
    }

	if (nPres != 0)
	    if (isAdded2[nPres-1]<0) isAdded2[nPres-1] = nClusters++;

    // Secondary clustering
    for (Int_t i=0; i<simLength; i+=2) { // use info from similarities2[] to merge preclusters
		for (Int_t m=0; m<nPres; m++)
			if (isAdded2[m] == similarities2[i])
				isAdded2[m] = similarities2[i+1];
    }
    for (Int_t i=0; i<nClusters; i++) {
		Int_t n = 0;
		for (Int_t j=0; j<nPres; j++)
			if (isAdded2[j]==i) n++;
		if (n > 0) {
			for (Int_t j=0; j<nPres; j++)
				if (isAdded2[j]==i) isAdded2[j] = clusterNr; // info in isAdded2: precluster "index" belongs to cluster "isAdded2[index]"
			clusterNr++;
		}
	}

	if (fVerbose>3){
		cout<<"Finished precluster assignment: "<<clusterNr<<" clusters identified\n"<<endl;
	}

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

	for (Int_t i=0; i<nPres; i++) {
		PndEmcCluster* thisPrecluster = (PndEmcCluster*) fPreclusterArray->At(i);
		if (isAdded2[i] < fClusterArray->GetEntriesFast() ) { // if isAdded2[i] is smaller than the current #clusters, the cluster to which this precluster belongs already exists, and we should add it to that
			PndEmcCluster* cluster= (PndEmcCluster*) fClusterArray->At(isAdded2[i]); // set pointer to the cluster we need, as indicated by isAdded2[]
			cluster->addCluster(thisPrecluster,fDigiArray); // add precluster i to the current cluster
		}
		else { // if not, make a new cluster
			PndEmcCluster* newCluster = new((*fClusterArray)[fClusterArray->GetEntriesFast()]) PndEmcCluster();
			newCluster->addCluster(thisPrecluster,fDigiArray);
		}
	}



// ----- FINISHED MERGING PRECLUSTERS ------------------------------------
	
	FinishClusters();

	CNtimer.Stop();
	Double_t CNrtime = CNtimer.RealTime();
	Double_t CNctime = CNtimer.CpuTime();
	DCtimer.Reset();

	CNtotRtime += CNrtime; // keep track how much time was spent on clustering
	CNtotCtime += CNctime;

	evtCounter++;
}


void PndEmcDistributedClustering::FinishPreclusters() {

	Int_t nCluster = fPreclusterArray->GetEntriesFast();
	if (fVerbose>3){
		cout<<nCluster<<")"<<endl;
	}
	for (Int_t i=0; i<nCluster; i++) {
		FinishPrecluster((PndEmcCluster*) (fPreclusterArray->At(i)));
	}	
}

void PndEmcDistributedClustering::FinishPrecluster(PndEmcCluster* precluster) {

	PndEmcClusterProperties clustProperties(*precluster, fDigiArray);

	TVector3 tmpprecpos = clustProperties.Where(fRecoPar->GetEmcClusterPosMethod(), fClusterPosParam);
	precluster->SetPosition(tmpprecpos);

	std::vector<Int_t> list = precluster->DigiList();

	Double_t Etot = 0;
	Double_t Emax = 0;
	Double_t maxDist = 0;
	Int_t Emax_idx = 0;

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

		Etot += thedigi->GetEnergy();
		if(thedigi->GetEnergy() > Emax) {
			Emax = thedigi->GetEnergy();
			Emax_idx = idx;
		}
		if (precluster->DistanceToCentre(thedigi)>maxDist) { // get cluster size (rough estimate)
			maxDist = precluster->DistanceToCentre(thedigi);
		}

	}
	precluster->SetRadius(maxDist); // new function added to PndEmcCluster (correction for crystal size already in function)
	precluster->SetEnergy(Etot);
	Double_t distanceToIP = tmpprecpos.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

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

	if (fVerbose>3) cout << "Assigned precluster properties: radius = " << precluster->radius() << ", energy = " << Etot << endl;

	// -------------- stop (end of precluster finding) -------------------------------------------------------
}

void PndEmcDistributedClustering::FinishClusters() {

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

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

void PndEmcDistributedClustering::FinishCluster(PndEmcCluster* cluster) {

	//----set energy and time etc. to the clusters
	
	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;
		}

	}

	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 PndEmcDistributedClustering::FinishTask() {
	if(fVerbose>2) {
		std::cout << "Processed in total " << digiCounter << " digis within " << evtCounter << " events." << std::endl;
	}
	cout << "[INFO\t] Preclustering (DC): real time: " << DCtotRtime << " s, cpu time: " << DCtotCtime << " s." << endl;
	cout << "[INFO\t] Clustering (CN): real time: " << CNtotRtime << " s, cpu time: " << CNtotCtime << " s." << endl;
}

void PndEmcDistributedClustering::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 PndEmcDistributedClustering::SetStorageOfData(Bool_t val) {
	fStoreClusters=val;
	return;
}



ClassImp(PndEmcDistributedClustering)