//-----------------------------------------------------------------------
// File and Version Information:
//      $Id: $
//
// Description:
//	Class EmcClusterLiloPos.
//
//	 Software developed for the BaBar Detector at the SLAC B-Factory.
// Adapted for the PANDA experiment at GSI		
//
// Author List:
//	Stephan Otto
// Dima Melnichuk - adaption for PANDA
//		
// Copyright Information:
//	Copyright (C) 2000	TU Dresden
//------------------------------------------------------------------------

#include "EmcClusterLiloPos.h"
#include "TwoCoordIndex.h"
#include "EmcCluster.h"
#include "EmcDigi.h"
#include "EmcXtal.h"		
#include "EmcStructure.h"		
#include "TVector3.h"

#include <iostream>
#include "assert.h"

EmcClusterLiloPos::EmcClusterLiloPos( const EmcCluster& toUse )
{
}

EmcClusterLiloPos::EmcClusterLiloPos( const EmcClusterLiloPos& other )
{
}

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

EmcClusterLiloPos::~EmcClusterLiloPos()
{
}

TVector3 
EmcClusterLiloPos::liloWhere( const EmcCluster *aClus )
{
	const double lOffsetParmA=3.6;//3.600; //4.0
	const double lOffsetParmB=1.594;
	const double lOffsetParmC=2.543;
	
	const double lClusEnergy=aClus->energy();
	
	assert(lClusEnergy!=0);
	
	const double lOffset=
		lOffsetParmA-lOffsetParmB*exp(-lOffsetParmC*lClusEnergy);
	
	TVector3 lLiloPoint( 1, 1, 1 );
	
	TVector3 lLinSum( 0, 0, 0 );
	TVector3 lLogSum( 0, 0, 0 );
	
	double lLogWeightSum=0;
	int lLogNum=0;
	
	bool lLogSecondTheta  = false;
	bool lLogSecondPhi    = false;
	int  lLogFirstTheta   = -666;
	int  lLogFirstPhi     = -666;
	
	EmcDigiPtrDict *lDigis = 
		(EmcDigiPtrDict *)aClus->itsMemberDigiMap();
	
	EmcDigiPtrDict::iterator lDigiIter = (*lDigis).begin();
	
	std::map<TwoCoordIndex*, EmcXtal*> tciXtalMap=EmcStructure::Instance("")->GetTciXtalMap();

	while(lDigiIter != (*lDigis).end())
	{
		EmcDigi *lDigi=lDigiIter->second;
		EmcXtal* xtal = tciXtalMap[lDigi->GetTCI()];
		const TVector3 tNormal=xtal->normalToFrontFace();
		TVector3 tmpPoint=lDigi->where();

		// TODO Consider forwars endcup separetly
		/*		if(lDigi->theta() > 3000)
		tmpPoint += *tNormal * 7.0 * 6.2;
		else*/
		tmpPoint += tNormal * 6.2;

		const TVector3 lDigiWhere=tmpPoint;

		const int lDigiTheta=lDigi->GetThetaInt();
		const int lDigiPhi=lDigi->GetPhiInt();
		const double lDigiEnergy=lDigi->GetEnergy();
		const double lLinWeight=lDigiEnergy/lClusEnergy;
		const double lLogWeight=lOffset+log(lLinWeight);
		
		lLinSum+=lLinWeight*lDigiWhere;
		
		if(lLogWeight>0)
		{
			lLogSum+=lLogWeight*lDigiWhere;
			lLogWeightSum+=lLogWeight;
			lLogNum++;
			
			if(lLogNum==1)
			{ 
					lLogFirstTheta=lDigiTheta;
					lLogFirstPhi=lDigiPhi;
			}
			else
			{
				if(!lLogSecondTheta&&lDigiTheta!=lLogFirstTheta) 
					lLogSecondTheta=true;
				if(!lLogSecondPhi&&lDigiPhi!=lLogFirstPhi) 
					lLogSecondPhi=true;
			}
		}
		++lDigiIter;
	}
  
	if(lLogNum>0) lLogSum*=1/lLogWeightSum;
	
	lLiloPoint.SetTheta(lLogSecondTheta?lLogSum.Theta():lLinSum.Theta());
	lLiloPoint.SetPhi(lLogSecondPhi?lLogSum.Phi():lLinSum.Phi());
	
	TwoCoordIndex *lTCI=EmcStructure::Instance("")->locateIndex(lLiloPoint.Theta(),lLiloPoint.Phi());
  
	assert(lTCI!=0);  
	
	const EmcXtal *lGeom=tciXtalMap[lTCI];
	
	// First, find out if the point is outside the crystal.
	
	const TVector3 lLiloVector = lLiloPoint;
	const TVector3 lCentre = lGeom->frontCentre();
	
	const TVector3 lNormal = lGeom->normalToFrontFace();
	const TVector3 lVector = lCentre - lLiloVector;
	const double lLength = lNormal.Dot(lVector);
	//const Hep3Vector lUnit = lLiloVector.unit();
	//const double lMag = (lNormal->dot(lCentre)/lNormal->dot(lUnit));  
	//cout<<" With original code, lilo-position "<<float(lMag)<<endl;

	if (lLength < 0.0)
	{
		// Point is outside crystal (not sure its correct for all cases?). 
		// Anyhow, project point back onto front-face.
		const TVector3 lUnit = lLiloVector.Unit();
		const double lMag = (lNormal.Dot(lCentre)/lNormal.Dot(lUnit));  
		lLiloPoint.SetMag(lMag);
	}
	else
	{
	// Don't project onto front-face. Keep the centroid inside the crytal.
	//cout<<" Point is inside crystal, keep as it is "<<endl;
		if (lLogNum > 1)
		{
			// Use logarithmic centroid position
			lLiloPoint.SetMag(lLogSum.Mag());
			//cout<<" With log centroid method, lilo-position "<<lLiloPoint.mag()<<endl;
		}
		else
		{
			// Use linear centroid position
			lLiloPoint.SetMag(lLinSum.Mag());
			//cout<<" Use lin centroid method, lilo-position "<<lLiloPoint.mag()<<endl;
		}
	
	}

	return lLiloPoint;
}

//ClassImp(EmcClusterLiloPos)