// $Id$
//
//    File: DEventRFBunch_factory_L3.cc
// Created: Thu Dec  3 17:27:55 EST 2009
// Creator: pmatt (on Linux ifarml6 2.6.18-128.el5 x86_64)
//

#include "DEventRFBunch_factory_L3.h"

using namespace jana;

//------------------
// init
//------------------
jerror_t DEventRFBunch_factory_L3::init(void)
{
	return NOERROR;
}

//------------------
// brun
//------------------
jerror_t DEventRFBunch_factory_L3::brun(jana::JEventLoop *locEventLoop, int32_t runnumber)
{
	DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
	DGeometry* locGeometry = locApplication->GetDGeometry(runnumber);

	dRFBunchFrequency = 4.008;

	double locTargetCenterZ;
	locGeometry->GetTargetZ(locTargetCenterZ);
	dTargetCenter.SetXYZ(0.0, 0.0, locTargetCenterZ);

	return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t DEventRFBunch_factory_L3::evnt(jana::JEventLoop *locEventLoop, uint64_t eventnumber)
{
	vector<const DRFTime*> locRFTimes;
	locEventLoop->Get(locRFTimes);
	if(locRFTimes.empty())
		return RESOURCE_UNAVAILABLE;

	double locRFHitTime = locRFTimes[0]->dTime;
	double locTimeVariance = locRFTimes[0]->dTimeVariance;

	//then figure out which RF bunch matches the most # tracks
		//need to try different sources of times: start with best quality
	int locBestRFBunchShift = 0;
	bool locNumParticleVotes = 0;
	vector<double> locTimes;

	if(Find_NeutralTimes(locEventLoop, locTimes)) //use neutral showers
		locNumParticleVotes = locTimes.size();
		locBestRFBunchShift = Find_BestRFBunchShift(locRFHitTime, locTimes);
	
	DEventRFBunch *locEventRFBunch = new DEventRFBunch;
	locEventRFBunch->dTime = locRFHitTime + dRFBunchFrequency*double(locBestRFBunchShift);
	locEventRFBunch->dTimeVariance = locTimeVariance;
	locEventRFBunch->dNumParticleVotes = locNumParticleVotes;
	_data.push_back(locEventRFBunch);

	return NOERROR;
}

bool DEventRFBunch_factory_L3::Find_NeutralTimes(JEventLoop* locEventLoop, vector<double>& locTimes)
{
	locTimes.clear();

	vector<const DBCALShower*> locBCALShowers;
	locEventLoop->Get(locBCALShowers);

	for(size_t loc_i = 0; loc_i < locBCALShowers.size(); ++loc_i)
	{
		if(locBCALShowers[loc_i]->E < 0.5) // require large energy shower
			continue;
		
		DVector3 locHitPoint(locBCALShowers[loc_i]->x, locBCALShowers[loc_i]->y, locBCALShowers[loc_i]->z);
		DVector3 locPath = locHitPoint - dTargetCenter;
		double locPathLength = locPath.Mag();
		if(!(locPathLength > 0.0))
			continue;

		double locFlightTime = locPathLength/29.9792458;
		double locHitTime = locBCALShowers[loc_i]->t;
		locTimes.push_back(locHitTime - locFlightTime);
	}

	vector<const DFCALShower*> locFCALShowers;
	locEventLoop->Get(locFCALShowers);

	for(size_t loc_i = 0; loc_i < locFCALShowers.size(); ++loc_i)
	{
		if(locFCALShowers[loc_i]->getEnergy() < 1.5) // require large energy shower
			continue;
		
		DVector3 locHitPoint = locFCALShowers[loc_i]->getPosition();
		DVector3 locPath = locHitPoint - dTargetCenter;
		double locPathLength = locPath.Mag();
		if(!(locPathLength > 0.0))
			continue;

		double locFlightTime = locPathLength/29.9792458;
		double locHitTime = locFCALShowers[loc_i]->getTime();
		locTimes.push_back(locHitTime - locFlightTime);
	}

	return (!locTimes.empty());
}

int DEventRFBunch_factory_L3::Find_BestRFBunchShift(double locRFHitTime, const vector<double>& locTimes)
{
	//then find the #beam buckets the RF time needs to shift to match it
	map<int, unsigned int> locNumRFBucketsShiftedMap; //key is # RF buckets the RF time is shifted to match the track, value is the # tracks at that shift
	int locBestRFBunchShift = 0;
	for(unsigned int loc_i = 0; loc_i < locTimes.size(); ++loc_i)
	{
		double locPropagatedTrackTime = locTimes[loc_i];
		//do manually: tricky to convert int to float...
		int locNumRFBucketsShifted = 0;
		double locTempRFHitTime = locRFHitTime;
		while((locTempRFHitTime - locPropagatedTrackTime) > (0.5*dRFBunchFrequency))
		{
			locTempRFHitTime -= dRFBunchFrequency;
			--locNumRFBucketsShifted;
		}
		while((locTempRFHitTime - locPropagatedTrackTime) < (-0.5*dRFBunchFrequency))
		{
			locTempRFHitTime += dRFBunchFrequency;
			++locNumRFBucketsShifted;
		}

		if(locNumRFBucketsShiftedMap.find(locNumRFBucketsShifted) == locNumRFBucketsShiftedMap.end())
			locNumRFBucketsShiftedMap[locNumRFBucketsShifted] = 1;
		else
			++(locNumRFBucketsShiftedMap[locNumRFBucketsShifted]);

		if(locNumRFBucketsShifted == locBestRFBunchShift)
			continue;

		unsigned int locBestNumTracks = locNumRFBucketsShiftedMap[locBestRFBunchShift];
		unsigned int locNumTracks = locNumRFBucketsShiftedMap[locNumRFBucketsShifted];
		if(locNumTracks > locBestNumTracks)
			locBestRFBunchShift = locNumRFBucketsShifted;
	}
	return locBestRFBunchShift;
}

//------------------
// erun
//------------------
jerror_t DEventRFBunch_factory_L3::erun(void)
{
	return NOERROR;
}

//------------------
// fini
//------------------
jerror_t DEventRFBunch_factory_L3::fini(void)
{
	return NOERROR;
}