//--------------------------------------------------------------------------
// File and Version Information:
// 	$Id:$
//
// Description:
//	Class PndEmcExpClusterSplitter.
//      Implementation of ClusterSplitter which splits
//      on the basis of exponential distance from the bump centroid.
//
// Environment:
//	Software developed for the BaBar Detector at the SLAC B-Factory.
//
// Adapted for the PANDA experiment at GSI
//
// Author List:
//      Phil Strother               
//
// Copyright Information:
//	Copyright (C) 1997               Imperial College
//
// Modified:
// M. Babai
//------------------------------------------------------------------------


//-----------------------
// This Class's Header --
//-----------------------
#include "PndEmcExpClusterSplitter.h"

//---------------
// C++ Headers --
//---------------
//#include <vector>
//#include <set>
//#include <map>
#include <iostream>

//-------------------------------
// Collaborating Class Headers --
//-------------------------------

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

#include "PndEmcClusterProperties.h"
#include "PndEmcXClMoments.h"

#include "PndEmcRecoPar.h"
#include "PndEmcGeoPar.h"
#include "PndEmcDigiPar.h"
#include "PndEmcStructure.h"
#include "PndEmcMapper.h"

#include "PndDetectorList.h"
#include "PndEmcTwoCoordIndex.h"
#include "PndEmcBump.h"
#include "PndEmcCluster.h"
#include "PndEmcDigi.h"
#include "PndEmcSharedDigi.h"
#include "PndEmcXtal.h"
#include "PndEmcDataTypes.h"
using std::endl;

//----------------
// Constructors --
//----------------

PndEmcExpClusterSplitter::PndEmcExpClusterSplitter(Int_t verbose):
fDigiArray(new TClonesArray()), fClusterArray(new TClonesArray()), fBumpArray(new TClonesArray()), fSharedDigiArray(new TClonesArray()), fGeoPar(new PndEmcGeoPar()), fDigiPar(new PndEmcDigiPar()), fRecoPar(new PndEmcRecoPar()), fPersistance(kTRUE), fMoliereRadius(0), fMoliereRadiusShashlyk(0), fExponentialConstant(0), fMaxIterations(0), fCentroidShift(0), fMaxBumps(0), fMinDigiEnergy(0), fClusterPosParam(), fVerbose(verbose)
{
  fClusterPosParam.clear();
}

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

PndEmcExpClusterSplitter::~PndEmcExpClusterSplitter()
{
// 	delete fGeoPar;
// 	delete fDigiPar;
// 	delete fRecoPar;
}

InitStatus PndEmcExpClusterSplitter::Init() {
  
  // Get RootManager
  FairRootManager* ioman = FairRootManager::Instance();
  if ( ! ioman ){
    cout << "-E- PndEmcMakeBump::Init: "
	 << "RootManager not instantiated!" << endl;
    return kFATAL;
  }
	
  // Geometry loading
  fGeoPar->InitEmcMapper();
  PndEmcStructure::Instance();
	
	// Get input array
	fDigiArray = (TClonesArray*) ioman->GetObject("EmcDigi");
	if ( ! fDigiArray ) {
		cout << "-W- PndEmcMakeCluster::Init: "
		<< "No PndEmcDigi array!" << endl;
		return kERROR;
	}
  
  fClusterArray = (TClonesArray*) ioman->GetObject("EmcCluster");
  if ( ! fClusterArray ) {
    cout << "-W- PndEmcMakeBump::Init: "
	 << "No PndEmcCluster array!" << endl;
    return kERROR;
  }
  
  fMoliereRadius=fRecoPar->GetMoliereRadius();// Mr in cm
  fMoliereRadiusShashlyk=fRecoPar->GetMoliereRadiusShashlyk();// Mr in cm
  fExponentialConstant=fRecoPar->GetExponentialConstant();
  // Energy fall off with distance from centre of cluster is
  // exp(-a*dist/Mr) Mr is the moliere radius.  dist is distance from
  // centre in cm, a is the above parameter.  The optimized value for
  // logarithimc cluster positioning is 2.5 For linear cluster
  // positioning a should be set to 1.5
  
  fMaxIterations=fRecoPar->GetMaxIterations();
  // Set the max number of iterations that the splitting algorithm is allowed 
  // to do before it decides that enough is enough.
  
  fCentroidShift=fRecoPar->GetCentroidShift();
  // Set the tolerance level to which it is said that a bump is 
  // said not to have moved since the last iteration.  
  // The default is a millimetre.
  fMaxBumps=fRecoPar->GetMaxBumps();
  // Set an upper limit on the number of bumps allowed in a cluster if it 
  // is to be split.
  
  fMinDigiEnergy=fRecoPar->GetMinDigiEnergy();
  // Set minimum SharedDigi energy to 20keV.
	
	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());
	}	

  // Create and register output array
  fBumpArray = new TClonesArray("PndEmcBump");
  ioman->Register("EmcBump","Emc",fBumpArray,fPersistance);
  
  fSharedDigiArray = new TClonesArray("PndEmcSharedDigi");
  ioman->Register("EmcSharedDigi","Emc",fSharedDigiArray,fPersistance);

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

}

void PndEmcExpClusterSplitter::Exec(Option_t* opt) 
{

	PndEmcMapper *fEmcMap=PndEmcMapper::Instance();
	// Reset output array
	if ( ! fBumpArray ) Fatal("Exec", "No Bump Array");
	fBumpArray->Delete();
	if ( ! fSharedDigiArray ) Fatal("Exec", "No Shared Digi Array");
	fSharedDigiArray->Delete();
	
	int nClusters = fClusterArray->GetEntriesFast();
	
	for (Int_t iCluster = 0; iCluster < nClusters; iCluster++){
		PndEmcCluster* theCluster = (PndEmcCluster*) fClusterArray->At(iCluster);

		int numberOfBumps = -1;
		numberOfBumps = (theCluster->LocalMaxMap()).size();
	
		if (numberOfBumps<=1 || numberOfBumps >= fMaxBumps){
			// Limit the max number of bumps in the cluster to 8 (default)
			// in this case, we clearly have a cluster, but no bumps to speak of.  
			// Make 1 bump with weights all equal to 1
		
			std::map<Int_t, Int_t>::const_iterator theDigiIterator;
			
			PndEmcBump* theNewBump = AddBump(); 
			theNewBump->MadeFrom(iCluster);
			theNewBump->SetLink(FairLink("EmcCluster", iCluster));

			for(theDigiIterator = theCluster->MemberDigiMap().begin();
			theDigiIterator != theCluster->MemberDigiMap().end(); ++theDigiIterator){
				PndEmcDigi *theDigi = (PndEmcDigi *) fDigiArray->At(theDigiIterator->second);
				PndEmcSharedDigi* sharedDigi=AddSharedDigi(theDigi, 1.0);
				Int_t iSharedDigi=fSharedDigiArray->GetEntriesFast()-1;
				theNewBump->addDigi(fSharedDigiArray,iSharedDigi);
			}
		}
		else {
			std::map<Int_t,Int_t> theMaximaDigis=theCluster->LocalMaxMap();
			std::map<Int_t,Int_t>::iterator theMaximaDigisIterator;
			std::map<PndEmcTwoCoordIndex*, TVector3*> theCentroidPoints;
			std::map<PndEmcTwoCoordIndex*, TVector3*> theMaximaPoints;
			std::map<PndEmcTwoCoordIndex*, PndEmcBump*> theIndexedBumps;
			std::map<PndEmcTwoCoordIndex*, TVector3*> theAllDigiPoints;
			
			std::map<Int_t, Int_t> theDigiDict = theCluster->MemberDigiMap();
			
			double totalEnergy=0;
			
			for(theMaximaDigisIterator = theMaximaDigis.begin(); 
			theMaximaDigisIterator != theMaximaDigis.end();
			++theMaximaDigisIterator){
				PndEmcDigi *theMaxDigi = (PndEmcDigi *) fDigiArray->At(theMaximaDigisIterator->second);
				
				Int_t detId = theMaximaDigisIterator->first;
				PndEmcTwoCoordIndex *theTCI = fEmcMap->GetTCI(detId);
				
				totalEnergy += theMaxDigi->GetEnergy();
				
				TVector3 *digiLocation = new TVector3(theMaxDigi->where());
				TVector3 *sameLocation = new TVector3(theMaxDigi->where());
				
				theMaximaPoints.insert(std::map<PndEmcTwoCoordIndex*, TVector3*>::value_type(theTCI, digiLocation));
				theCentroidPoints.insert(std::map<PndEmcTwoCoordIndex*, TVector3*>::value_type(theTCI, sameLocation));
			}
	
			std::map<Int_t,Int_t>::iterator theAllDigisIterator;		
			// This loop works out the location of all the digis in the cluster 
		
			for( theAllDigisIterator = theDigiDict.begin(); theAllDigisIterator != theDigiDict.end();
			++theAllDigisIterator){
				Int_t detId = theAllDigisIterator->first;
				PndEmcTwoCoordIndex *theTCI = fEmcMap->GetTCI(detId);
				PndEmcDigi *theDigi = (PndEmcDigi *) fDigiArray->At(theAllDigisIterator->second);
				TVector3 *digiLocation = new TVector3(theDigi->where());
				theAllDigiPoints.insert(std::map<PndEmcTwoCoordIndex*, TVector3*>::value_type( theTCI, digiLocation ));
			}
	
			theMaximaDigisIterator = theMaximaDigis.begin();
			
			// Now we can create the EmcBumps
			
			// The algorithm is as follows: We will index each bump by its
			// maximum digi's PndEmcTwoCoordIndex.  We will set up a list of
			// bump centroids which to start with will be synonymous with the
			// location of the maxima.  We then apportion a weight to each
			// digi, according to its distance from the centroids.  We then
			// construct the bumps according to these weights, which will
			// presumably give a different set of centroids.  This is repeated
			// until the centroids are static within tolerance, or we reach
			// the maximum number of iterations.
			
			Int_t iterations = 0;
			
			Double_t averageCentroidShift;
		
		do {
			if (fVerbose>=3){
				std::cout<<"iteration No "<<iterations<<std::endl;
			}
			averageCentroidShift=0.0;
			
	      // First clean up the old bumps
			std::map<PndEmcTwoCoordIndex*, PndEmcBump*>::iterator theBumpKiller = theIndexedBumps.begin();
	      while(theBumpKiller != theIndexedBumps.end()){
	        PndEmcBump *theBump = theBumpKiller->second;
	        delete theBump;
	        ++theBumpKiller;
	      }
	      theIndexedBumps.clear();
	      			
			// Then loop over all the maxima and assign weights accordingly
			for (theMaximaDigisIterator = theMaximaDigis.begin();
			theMaximaDigisIterator != theMaximaDigis.end();
			++theMaximaDigisIterator) {
				Int_t detId = theMaximaDigisIterator->first;
				PndEmcTwoCoordIndex *theCurrentMaximaTCI = fEmcMap->GetTCI(detId);
			
				if (fVerbose>=3){
				std::cout<<"***************** current maximum: theta = "<<theCurrentMaximaTCI->XCoord()
						<<", phi = "<<theCurrentMaximaTCI->YCoord()<<"*********"<<std::endl;
				}
			
				// Create the bump which will correspond to this digi maxima
				PndEmcBump* theNewBump = new PndEmcBump();
				theNewBump->MadeFrom(iCluster);
				theIndexedBumps.insert(std::map<PndEmcTwoCoordIndex*,
					PndEmcBump*>::value_type(theCurrentMaximaTCI, theNewBump));

			
				// Now we will look over all the digis and add each of them
				// to this Bump with an appropriate weight
			
				for (theAllDigisIterator = theDigiDict.begin();
					theAllDigisIterator != theDigiDict.end();++theAllDigisIterator) {
					PndEmcDigi *theCurrentDigi = (PndEmcDigi *) fDigiArray->At(theAllDigisIterator->second);
					PndEmcTwoCoordIndex *theCurrentTCI = theCurrentDigi->GetTCI();
			
					Double_t weight;
			
				// We are on the first pass and the digi is not a local max, or we are not on the
				// first pass.   Assign a weight according to the distance from the centroid position.
			
					Double_t myEnergy = 0;
					Double_t myDistance = 0;
					
					// Now share the digi out according to its distance from the maxima, 
					// and the maxima energies
					
					Double_t totalDistanceEnergy=0;
		
					//Moliere Radius for Shashlyk is different
					Double_t MoliereRadius;
					if(theCurrentDigi->GetModule() == 5)
						MoliereRadius = fMoliereRadiusShashlyk;
					else
						MoliereRadius = fMoliereRadius;

					std::map<PndEmcTwoCoordIndex*,TVector3*>::iterator theMaxPointsIterator;
					
					for(theMaxPointsIterator = theCentroidPoints.begin(); theMaxPointsIterator != theCentroidPoints.end();++theMaxPointsIterator){
						PndEmcTwoCoordIndex *theMaxPointsTCI =theMaxPointsIterator->first; 
				
						TVector3 *theMaxPoint = theMaxPointsIterator->second;
						TVector3 *theCurrentDigiPoint = theAllDigiPoints.find(theCurrentTCI)->second;
						
						Double_t theDistance;
						
						// This next bit just checks to see if the maxima point in
						// hand is the same as the crystal from which we are
						// trying to find distance - just an FP trap really.
						
						if ((*theCurrentTCI)==(*theMaxPointsTCI)){
							theDistance=0.0;
						} else {
							TVector3 distance( *theMaxPoint - *theCurrentDigiPoint);
							
							theDistance = distance.Mag();
						}
						
						if (*theCurrentMaximaTCI == *(theMaxPointsTCI)){
						// i.e. the maximum we are trying to find the distance from is 
						// the one for which we are currently trying to make a bump
						myDistance = theDistance;
						Int_t iCurentMaxDigi = (theDigiDict.find(theMaxPointsTCI->Index()))->second;
						myEnergy = ((PndEmcDigi *) fDigiArray->At(iCurentMaxDigi))->GetEnergy();
						}
				
						Int_t iMaxPoint = (theDigiDict.find(theMaxPointsTCI->Index()))->second;
						totalDistanceEnergy += ((PndEmcDigi *) fDigiArray->At(iMaxPoint))->GetEnergy() *
							exp(-fExponentialConstant * theDistance/MoliereRadius);
					}
		
					if(totalDistanceEnergy > 0.0)
						weight = myEnergy*exp(-fExponentialConstant* 
								myDistance/MoliereRadius) / ( totalDistanceEnergy);
					else 
						weight=0;
					
					if (fVerbose>=3){
						std::cout<<"\t digi theta = "<<theCurrentDigi->GetTCI()->XCoord()
							<<", phi = "<<theCurrentDigi->GetTCI()->YCoord()<<std::endl;
						std::cout<<"energy = "<<theCurrentDigi->GetEnergy()<<", weight = "<< weight<<std::endl;
					}
					PndEmcSharedDigi* sharedDigi= AddSharedDigi( theCurrentDigi, weight );
				
					if (fVerbose>=3){
						std::cout<<"shared digi energy = "<<sharedDigi->GetEnergy()<<std::endl;
					}
					
					Int_t iSharedDigi=fSharedDigiArray->GetEntriesFast()-1;
					if( sharedDigi->GetEnergy() > fMinDigiEnergy){
						theNewBump->addDigi( fSharedDigiArray, iSharedDigi );
					} else {
						fSharedDigiArray->RemoveAt(iSharedDigi);
						fSharedDigiArray->Compress();
					}
				}
		
				// Compute the shift of the centroid we have just calculated
				TVector3 *theOldCentroid = theCentroidPoints.find(theCurrentMaximaTCI)->second;
				
				PndEmcClusterProperties clusterProperties(*theNewBump, fSharedDigiArray);
				
				TVector3 newbumppos =  clusterProperties.Where(fRecoPar->GetEmcClusterPosMethod(), fClusterPosParam);
				theNewBump->SetPosition(newbumppos);
				TVector3 centroidShift(*theOldCentroid - newbumppos);
				averageCentroidShift+=centroidShift.Mag();
			}
			
			averageCentroidShift/=(Double_t)numberOfBumps;
			
			// Put the new centroids in the list of centroid points,
			// remembering to delete the old ones.
			std::map<PndEmcTwoCoordIndex*, TVector3*>::iterator 
			theCentroidPointsIterator = theCentroidPoints.begin();
			for(theCentroidPointsIterator = theCentroidPoints.begin();
			theCentroidPointsIterator != theCentroidPoints.end(); 
			++theCentroidPointsIterator) {
				delete theCentroidPointsIterator->second;
			}
			theCentroidPoints.clear();
     
			std::map<PndEmcTwoCoordIndex*, PndEmcBump*>::iterator theIndexedBumpsIterator;
	      for(theIndexedBumpsIterator = theIndexedBumps.begin();
				theIndexedBumpsIterator != theIndexedBumps.end(); ++theIndexedBumpsIterator){
	        TVector3 *theNewCentroid = new TVector3((theIndexedBumpsIterator->second)->where());
	        theCentroidPoints.insert(std::map<PndEmcTwoCoordIndex*, TVector3*>::value_type(theIndexedBumpsIterator->first,theNewCentroid));
			}

			iterations++;
			
		} while (iterations < fMaxIterations && averageCentroidShift > fCentroidShift);
		//End of do loop
		
		// Finally append the new bumps to the TClonesArray.
		std::map<PndEmcTwoCoordIndex*, PndEmcBump*>::iterator theBumpsIterator;
		PndEmcBump *theBump;
		
		for(theBumpsIterator = theIndexedBumps.begin(); theBumpsIterator != theIndexedBumps.end();++theBumpsIterator){
			theBump = theBumpsIterator->second;
			Int_t size_ba = fBumpArray->GetEntriesFast();
			PndEmcBump* theNextBump = new((*fBumpArray)[size_ba]) PndEmcBump(*(theBump));
			if (fVerbose>0)
				std::cout << "Bump Created!" << std::endl;
			theNextBump->SetLink(FairLink("EmcCluster", iCluster));
		}
		
		std::map<PndEmcTwoCoordIndex*,TVector3*>::iterator theGrimReaper = theMaximaPoints.begin();
		for(theGrimReaper = theMaximaPoints.begin();
		theGrimReaper != theMaximaPoints.end();
		++theGrimReaper){
			delete theGrimReaper->second;
		}
		theMaximaPoints.clear();
		
		for(theGrimReaper  = theAllDigiPoints.begin();
		theGrimReaper != theAllDigiPoints.end();
		++theGrimReaper){
			delete theGrimReaper->second;
		}
		theAllDigiPoints.clear();
		
		for(theGrimReaper = theCentroidPoints.begin(); 
		theGrimReaper != theCentroidPoints.end();
		++theGrimReaper){
			delete theGrimReaper->second;
		}
		theCentroidPoints.clear();
	}
	
	Int_t nBumps = (theCluster->LocalMaxMap()).size();
	theCluster->SetNBumps(nBumps);
	}
  
	// At that moment internal state fEnergy and fWhere of Clusters are
	// not initialized, the following make it possible to see energy and
	// position from output root file
	Int_t nBump = fBumpArray->GetEntriesFast();

	for (Int_t i=0; i<nBump; i++){
		PndEmcBump *tmpbump = (PndEmcBump*) fBumpArray->At(i);
		PndEmcClusterProperties clusterProperties(*tmpbump, fSharedDigiArray);
		
		if (!tmpbump->IsEnergyValid())
			tmpbump->SetEnergy(clusterProperties.Energy());
		if (!tmpbump->IsPositionValid())
			tmpbump->SetPosition(clusterProperties.Where(fRecoPar->GetEmcClusterPosMethod(), fClusterPosParam));
		PndEmcXClMoments xClMoments(*tmpbump, fSharedDigiArray);
		tmpbump->SetZ20(xClMoments.AbsZernikeMoment(2, 0, 15));
		tmpbump->SetZ53(xClMoments.AbsZernikeMoment(5, 3, 15));
		tmpbump->SetLatMom(xClMoments.Lat());

	}
  
  if (fVerbose>=1){
    std::cout<<"PndEmcExpClusterSplitter:: Number of clusters = "<<nClusters<<std::endl;
    std::cout<<"PndEmcExpClusterSplitter:: Number of bumps = "<<nBump<<std::endl;
  }

}

PndEmcBump* PndEmcExpClusterSplitter::AddBump(){
  TClonesArray& clref = *fBumpArray;
  Int_t size = clref.GetEntriesFast();
  return new(clref[size]) PndEmcBump();
}

PndEmcSharedDigi* PndEmcExpClusterSplitter::AddSharedDigi(PndEmcDigi* digi, Double_t weight){
  TClonesArray& clref = *fSharedDigiArray;
  Int_t size = clref.GetEntriesFast();
  return new(clref[size]) PndEmcSharedDigi(*digi, weight);
}

void PndEmcExpClusterSplitter::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 digitisation 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");
}

ClassImp(PndEmcExpClusterSplitter)