/* 
 * File:   TBCluster.cc
 * Author: daniel
 * 
 * Created on 9. Februar 2014, 14:02
 */

#include "TBEvent.h"

// TBCluster
TBCluster::TBCluster()
{
	TBCluster::fMaskedCluster = kUnknown;
	TBCluster::fEdgePixelCluster = kUnknown;
	TBCluster::fCentralPixelCluster = kUnknown;
	TBCluster::fMaskedPixelCluster = kUnknown;
	TBCluster::fClusterGood = kUnknown;
}

TBCluster::~TBCluster()
{
	TBCluster::hits.clear();
}

// static
int TBCluster::getSumToT(const TBCluster* tbcluster)
{
	if(tbcluster == NULL)
	{
		return 0;
	}
	
	int sumToT = 0;
	for(auto hit: tbcluster->hits)
	{
		sumToT += hit->tot;
	}
	
	return sumToT;
}

int TBCluster::getSumCharge(const TBCluster* tbcluster, const TBEvent* event)
{
	double sumCharge = 0.0;
	for(auto hit: tbcluster->hits)
	{
		sumCharge += (!event->dut->totcalib->hasCalib()) ? hit->tot : event->dut->q(hit->tot, hit->col, hit->row);
	}

	return (int) sumCharge;
}

TBCluster* TBCluster::getMaxToTCluster(const std::vector<TBCluster*>* clusters)
{
	int maxToT = 0;
	TBCluster* maxToTCluster = NULL;
	for(auto tbcluster: (*clusters))
	{
		int tmpToT = TBCluster::getSumToT(tbcluster);
		if(tmpToT > maxToT)
		{
			maxToT = tmpToT;
			maxToTCluster = tbcluster;
		}
	}
	
	return maxToTCluster;
}

TBHit* TBCluster::getMaxToTCell(const TBCluster* tbcluster)
{
	int maxToT = 0;
	TBHit* maxToTHit = NULL;
	for(auto hit: tbcluster->hits)
	{
		if(hit->tot < maxToT)
		{ 
			continue;
		}
		maxToT = hit->tot;
		maxToTHit = hit;
	}
	return maxToTHit;
}

TBCluster* TBCluster::getMatchedCluster(const std::vector<TBCluster*>* clusters, const TBTrack* tbtrack, const TBEvent* event, double matchX /* =  0.0 */, double matchY /* = 0.0 */)
{
	TBCluster* matchCluster = NULL;
	double x, y;
	double uncertaintyX, uncertaintyY;

	// todo get match radius, nach x & y getrennt oder mit einem Radius?
	// ueberarbeiten !!!!!!!!!
	if(matchX == 0.0)
	{
		matchX = std::fabs(event->dut->getMatchX());
	}
	
	if(matchY == 0.0)
	{
		matchY = std::fabs(event->dut->getMatchY());
	}
	
	double minDis = std::numeric_limits<double>::max();
	
	for(auto cluster: *clusters)
	{
		TBCluster::getClusterCenter(cluster, event, &x, &y, &uncertaintyX, &uncertaintyY);
		double disX = std::fabs(tbtrack->trackX - x);
		double disY = std::fabs(tbtrack->trackY - y);
		if(disX < uncertaintyX + matchX and disY < uncertaintyY + matchY)
		{
			double dis = std::sqrt(disX*disX + disY*disY);
			if(dis < minDis)
			{
				minDis = dis;
				matchCluster = cluster;
			}
		}
	}
	return matchCluster;
}

bool TBCluster::isMatch(const TBCluster* tbcluster, const TBTrack* tbtrack, const TBEvent* event)
{
	double x, y;
	double uncertaintyX, uncertaintyY;
	
	// todo get match radius, nach x & y getrennt oder mit einem Radius?
	// ueberarbeiten !!!!!!!!!
	 double matchX = std::fabs(event->dut->getMatchX());
	 double matchY = std::fabs(event->dut->getMatchY());
	
	TBCluster::getClusterCenter(tbcluster, event, &x, &y, &uncertaintyX, &uncertaintyY);
	if(std::fabs(tbtrack->trackX - x) < uncertaintyX + matchX)
	{
		if (std::fabs(tbtrack->trackY - y) < uncertaintyY + matchY)
		{
			return true;
		}
	}
	
	return false;
}

double TBCluster::getClusterCenter(const TBCluster* tbcluster, const TBEvent* event, double* x, double* y, double* uncertaintyX, double* uncertaintyY)
{
	*x = TBCluster::getChargeWeightedX(tbcluster, event);
	*y = TBCluster::getChargeWeightedY(tbcluster, event);
	// todo add uncertainty
	*uncertaintyX = 0.0;
	*uncertaintyY = 0.0;

	//todo (Histogrammmethode): Histogramme der jeweiligen Pixel nach passenden Positionen (dut-Koordinatensystem) durchsuchen
	// Ueberschneidungen der verschiedenen Pixels feststellen?! --> gewichten?
	// in Fehlerberechnung & Rückgabe
	return 0.0;
}

/*bool cluster::isCentralRegion(const event::cluster_t &cluster, const Event &event) {  
  // check, if every hit of the cluster is central
  for(int ii = 0; ii < cluster.size(); ii++) {
    if (tbutils::inCentralRegion(cluster[ii], event) == false){
      return false;
    }
  }
  return true;
}*/

double TBCluster::getUnWeightedX(const TBCluster* tbcluster, const TBEvent* event)
{
	double pos = 0.0;
	for (auto hit: tbcluster->hits)
	{
		pos += event->dut->pixelArray[hit->col][hit->row].getGeometricCenterX();
	}
	return (pos/ ((double) tbcluster->hits.size()));
}

double TBCluster::getUnWeightedY(const TBCluster* tbcluster, const TBEvent* event)
{
	double pos = 0.0;
	for (auto hit: tbcluster->hits)
	{
		pos += event->dut->pixelArray[hit->col][hit->row].getGeometricCenterY();
	}
	return (pos/ ((double) tbcluster->hits.size()));
}

double TBCluster::getChargeWeightedX(const TBCluster* tbcluster, const TBEvent* event)
{
	double w = 0.0;
	double pos = 0.0;
	for(auto hit: tbcluster->hits)
	{
		 int col = hit->col;
		 int row = hit->row;
		 int tot = hit->tot;
		 double weight = (!event->dut->totcalib->hasCalib() || event->dut->q(tot, col, row) == 0) ? (double) tot : event->dut->q(tot, col, row);
		 w += weight;
		 pos += event->dut->pixelArray[col][row].getGeometricCenterX() * weight;
	}
	if(w == 0)
	{
		std::cerr << "[ TBCluster ]: Error in getChargeWeightedX(): Cluster with 0 tot detected. You have to check that this does not occur before calling getChargeWeightedX()." << std::endl;
		exit(-1);
	}
	return (pos/w);
}

double TBCluster::getChargeWeightedY(const TBCluster* tbcluster, const TBEvent* event)
{
	double w = 0.0;
	double pos = 0.0;
	for(auto hit: tbcluster->hits)
	{
		 int col = hit->col;
		 int row = hit->row;
		 int tot = hit->tot;
		 double weight = (!event->dut->totcalib->hasCalib() || event->dut->q(tot, col, row) == 0) ? (double) tot : event->dut->q(tot, col, row);
		 w += weight;
		 pos += event->dut->pixelArray[col][row].getGeometricCenterY() * weight;
	}
	if(w == 0)
	{
		std::cerr << "[ TBCluster ]: Error in getChargeWeightedY(): Cluster with 0 tot detected. You have to check that this does not occur before calling getChargeWeightedY()." << std::endl;
		exit(-1);
	}
	return (pos/w);
}

double TBCluster::getClusterCenterColByDigitalHeadTail1D(const TBCluster* tbcluster, const TBEvent* event)
{
	int head = event->dut->getNcols()-1;
	int tail = 0;
	
	for(auto hit: tbcluster->hits)
	{
		// find head
		if(hit->col < head)
		{
			head = hit->col;
		}
		// find tail
		if(hit->col > tail)
		{
			tail = hit->col;
		}
	}
	
	return ((head+tail)/2);
}

double TBCluster::getClusterCenterXByDigitalHeadTail1D(const TBCluster* tbcluster, const TBEvent* event)
{
	DUT* dut = event->dut;
	double centerCol = TBCluster::getClusterCenterColByDigitalHeadTail1D(tbcluster, event);
	double centerRow = TBCluster::getClusterCenterRowByDigitalHeadTail1D(tbcluster, event);
	
	double centerX = (dut->pixelArray[(int) std::floor(centerCol)][(int) std::floor(centerRow)].getGeometricCenterX() + dut->pixelArray[(int) std::ceil(centerCol)][(int) std::ceil(centerRow)].getGeometricCenterX())/2.0;
	
	return centerX;
}

double TBCluster::getClusterCenterRowByDigitalHeadTail1D(const TBCluster* tbcluster, const TBEvent* event)
{
	int head = event->dut->getNrows()-1;
	int tail = 0;
	
	for(auto hit: tbcluster->hits)
	{
		// find head
		if(hit->row < head)
		{
			head = hit->row;
		}
		// find tail
		if(hit->row > tail)
		{
			tail = hit->row;
		}
	}
	
	return ((head+tail)/2);
}

double TBCluster::getClusterCenterYByDigitalHeadTail1D(const TBCluster* tbcluster, const TBEvent* event)
{
	DUT* dut = event->dut;
	double centerCol = TBCluster::getClusterCenterColByDigitalHeadTail1D(tbcluster, event);
	double centerRow = TBCluster::getClusterCenterRowByDigitalHeadTail1D(tbcluster, event);
	
	double centerY = (dut->pixelArray[(int) std::floor(centerCol)][(int) std::floor(centerRow)].getGeometricCenterY() + dut->pixelArray[(int) std::ceil(centerCol)][(int) std::ceil(centerRow)].getGeometricCenterY())/2.0;
	
	return centerY;
}

double TBCluster::getClusterCenterColByAnalogHeadTail1D(const TBCluster* tbcluster, const TBEvent* event, int signal_0, double t, double theta)
{
	// signal_0: most probable signal released by a particle at 0 degree
	// t: 'detector thickness' or rather active thickness
	// theta: inclination
	// Beware that several presumptions (homogeneous charge collection in active zone, active zone thickness, inclination angle...) are made that are not necessarily true.
	DUT* dut = event->dut;
	double result;
	
	int head = event->dut->getNcols()-1;
	int tail = 0;
	
	int signal_head_pixel;
	int signal_tail_pixel;
	
	double midPixelPitchX = 0.0;
	
	for(auto hit: tbcluster->hits)
	{
		// find head
		if(hit->col < head)
		{
			head = hit->col;
			signal_head_pixel = hit->tot;
		}
		// find tail
		if(hit->col > tail)
		{
			tail = hit->col;
			signal_tail_pixel = hit->tot;
		}
		
		midPixelPitchX += dut->pixelArray[hit->col][hit->row].getGeometry()->getMaxPitchX();
	}
	
	midPixelPitchX /= (double) tbcluster->hits.size();
	double signal_theta = (signal_0 * midPixelPitchX)/(t * std::sin(theta));

	result = TBCluster::getClusterCenterColByDigitalHeadTail1D(tbcluster, event) + midPixelPitchX * (std::min((double) signal_head_pixel, signal_theta)-std::min((double) signal_tail_pixel, signal_theta))/(2*signal_theta);

	return result;
}

double TBCluster::getClusterCenterXByAnalogHeadTail1D(const TBCluster* tbcluster, const TBEvent* event, int signal_0, double t, double theta, double phi)
{
	DUT* dut = event->dut;
	double centerCol = TBCluster::getClusterCenterColByAnalogHeadTail1D(tbcluster, event, signal_0, t, theta);
	double centerRow = TBCluster::getClusterCenterRowByAnalogHeadTail1D(tbcluster, event, signal_0, t, phi);
	
	double centerX = (dut->pixelArray[(int) std::floor(centerCol)][(int) std::floor(centerRow)].getGeometricCenterX() + dut->pixelArray[(int) std::ceil(centerCol)][(int) std::ceil(centerRow)].getGeometricCenterX())/2.0;
	
	return centerX;
}

double TBCluster::getClusterCenterRowByAnalogHeadTail1D(const TBCluster* tbcluster, const TBEvent* event, int signal_0, double t, double phi)
{
	// signal_0: most probable signal released by a particle at 0 degree
	// t: 'detector thickness' or rather active thickness
	// phi: inclination
	// Beware that several presumptions (homogeneous charge collection in active zone, active zone thickness, inclination angle...) are made that are not necessarily true.
	DUT* dut = event->dut;
	double result;
	
	int head = event->dut->getNrows()-1;
	int tail = 0;
	
	int signal_head_pixel;
	int signal_tail_pixel;
	
	double midPixelPitchY = 0.0;
	
	for(auto hit: tbcluster->hits)
	{
		// find head
		if(hit->row < head)
		{
			head = hit->row;
			signal_head_pixel = hit->tot;
		}
		// find tail
		if(hit->row > tail)
		{
			tail = hit->row;
			signal_tail_pixel = hit->tot;
		}
		
		midPixelPitchY += dut->pixelArray[hit->col][hit->row].getGeometry()->getMaxPitchY();
	}
	
	midPixelPitchY /= (double) tbcluster->hits.size();
	double signal_phi = (signal_0 * midPixelPitchY)/(t * std::sin(phi));

	result = TBCluster::getClusterCenterRowByDigitalHeadTail1D(tbcluster, event) + midPixelPitchY * (std::min((double) signal_head_pixel, signal_phi) - std::min((double) signal_tail_pixel, signal_phi))/(2*signal_phi);

	return result;
}

double TBCluster::getClusterCenterYByAnalogHeadTail1D(const TBCluster* tbcluster, const TBEvent* event, int signal_0, double t, double theta, double phi)
{
	DUT* dut = event->dut;
	double centerCol = TBCluster::getClusterCenterColByAnalogHeadTail1D(tbcluster, event, signal_0, t, theta);
	double centerRow = TBCluster::getClusterCenterRowByAnalogHeadTail1D(tbcluster, event, signal_0, t, phi);
	
	double centerY = (dut->pixelArray[(int) std::floor(centerCol)][(int) std::floor(centerRow)].getGeometricCenterY() + dut->pixelArray[(int) std::ceil(centerCol)][(int) std::ceil(centerRow)].getGeometricCenterY())/2.0;
	
	return centerY;
}

double TBCluster::getClusterCenterColByOneSidedAnalogHeadTail1D(const TBCluster* tbcluster, const TBEvent* event, int headSide, int signal_0, double t, double theta)
{
	// headSide: -1 for left and 1 for right (default -1)
	// signal_0: most probable signal released by a particle at 0 degree
	// t: 'detector thickness' or rather active thickness
	// theta: inclination
	// Beware that several presumptions (homogeneous charge collection in active zone, active zone thickness, inclination angle...) are made that are not necessarily true.
	DUT* dut = event->dut;
	double result;

	int signal_head_pixel;
	int head;
	
	double midPixelPitchX = 0.0;
	
	switch(headSide)
	{
		case 1:
			head = 0;
			break;
		default:
		case -1:
			headSide = -1;
			head = dut->getNcols()-1;
			break;
	}

	//find head
	// FIXME: Determine which side is supposed to be the head side.
	for(auto hit: tbcluster->hits)
	{
		if(headSide == -1 and hit->col < head)
		{
			head = hit->col;
			signal_head_pixel = hit->tot;
		}
		if(headSide == 1 and hit->col > head)
		{
			head = hit->col;
			signal_head_pixel = hit->tot;
		}
		
		midPixelPitchX += dut->pixelArray[hit->col][hit->row].getGeometry()->getMaxPitchX();
	}
	midPixelPitchX /= tbcluster->hits.size();
	
	double signal_theta = (signal_0 * midPixelPitchX)/(t * std::sin(theta));

	result = head + 0.5 * midPixelPitchX - midPixelPitchX * signal_head_pixel/signal_theta + 0.5 * t * std::tan(theta);

	return result;
}

double TBCluster::getClusterCenterXByOneSidedAnalogHeadTail1D(const TBCluster* tbcluster, const TBEvent* event, int headSideCol, int headSideRow, int signal_0, double t, double theta, double phi)
{
	if(headSideCol != -1 or headSideCol != 1)
	{
		headSideCol = -1;
	}
	if(headSideRow != -1 or headSideRow != 1)
	{
		headSideRow = -1;
	}
	
	DUT* dut = event->dut;
	double centerCol = TBCluster::getClusterCenterColByOneSidedAnalogHeadTail1D(tbcluster, event, headSideCol, signal_0, t, theta);
	double centerRow = TBCluster::getClusterCenterRowByOneSidedAnalogHeadTail1D(tbcluster, event, headSideRow, signal_0, t, phi);
	
	double centerX = (dut->pixelArray[(int) std::floor(centerCol)][(int) std::floor(centerRow)].getGeometricCenterX() + dut->pixelArray[(int) std::ceil(centerCol)][(int) std::ceil(centerRow)].getGeometricCenterX())/2.0;
	
	return centerX;
}

double TBCluster::getClusterCenterRowByOneSidedAnalogHeadTail1D(const TBCluster* tbcluster, const TBEvent* event, int headSide, int signal_0, double t, double phi)
{
	// headSide: -1 for top and 1 for bottom (default -1)
	// signal_0: most probable signal released by a particle at 0 degree
	// t: 'detector thickness' or rather active thickness
	// phi: inclination
	// Beware that several presumptions (homogeneous charge collection in active zone, active zone thickness, inclination angle...) are made that are not necessarily true.
	DUT* dut = event->dut;
	double result;

	int signal_head_pixel;
	int head;
	
	double midPixelPitchY = 0.0;
	
	switch(headSide)
	{
		case 1:
			head = 0;
			break;
		default:
		case -1:
			headSide = -1;
			head = dut->getNrows()-1;
			break;
	}

	//find head
	// FIXME: Determine which side is supposed to be the head side.
	for(auto hit: tbcluster->hits)
	{
		if(headSide == -1 and hit->row < head)
		{
			head = hit->row;
			signal_head_pixel = hit->tot;
		}
		if(headSide == 1 and hit->row > head)
		{
			head = hit->row;
			signal_head_pixel = hit->tot;
		}
		
		midPixelPitchY += dut->pixelArray[hit->col][hit->row].getGeometry()->getMaxPitchY();
	}
	midPixelPitchY /= tbcluster->hits.size();
	
	double signal_phi = (signal_0 * midPixelPitchY)/(t * std::sin(phi));

	result = head + 0.5 * midPixelPitchY - midPixelPitchY * signal_head_pixel/signal_phi + 0.5 * t * std::tan(phi);

	return result;
}

double TBCluster::getClusterCenterYByOneSidedAnalogHeadTail1D(const TBCluster* tbcluster, const TBEvent* event, int headSideCol, int headSideRow, int signal_0, double t, double theta, double phi)
{
	if(headSideCol != -1 or headSideCol != 1)
	{
		headSideCol = -1;
	}
	if(headSideRow != -1 or headSideRow != 1)
	{
		headSideRow = -1;
	}
	
	DUT* dut = event->dut;
	double centerCol = TBCluster::getClusterCenterColByOneSidedAnalogHeadTail1D(tbcluster, event, headSideCol, signal_0, t, theta);
	double centerRow = TBCluster::getClusterCenterRowByOneSidedAnalogHeadTail1D(tbcluster, event, headSideRow, signal_0, t, phi);
	
	double centerY = (dut->pixelArray[(int) std::floor(centerCol)][(int) std::floor(centerRow)].getGeometricCenterY() + dut->pixelArray[(int) std::ceil(centerCol)][(int) std::ceil(centerRow)].getGeometricCenterY())/2.0;
	
	return centerY;
}

double TBCluster::getInClusterMaxPitchX(const TBCluster* tbcluster, const TBEvent* event)
{
	DUT* dut = event->dut;
	double maxPixelPitchX = 0;
	double pitchX;
	for(auto hit: tbcluster->hits)
	{
		pitchX = dut->pixelArray[hit->col][hit->row].getGeometry()->getMaxPitchX();
		if(pitchX > maxPixelPitchX)
		{
			maxPixelPitchX = pitchX;
		}
	}
	return maxPixelPitchX;
}
  
double TBCluster::getInClusterMaxPitchY(const TBCluster* tbcluster, const TBEvent* event)
{
	DUT* dut = event->dut;
	double maxPixelPitchY = 0;
	double pitchY;
	for(auto hit: tbcluster->hits)
	{
		pitchY = dut->pixelArray[hit->col][hit->row].getGeometry()->getMaxPitchY();
		if(pitchY > maxPixelPitchY)
		{
			maxPixelPitchY = pitchY;
		}
	}
	return maxPixelPitchY;
}

int TBCluster::spannedRows(const TBCluster* tbcluster, const TBEvent* event)
{
	DUT* dut = event->dut;
	int min = dut->getNrows()-1;
	int max = 0;
	for(auto hit: tbcluster->hits)
	{
		if(hit->row < min)
		{
			min = hit->row;
		}
		if(hit->row > max)
		{
			max = hit->row;
		}
	}
	return max - min + 1;
}

int TBCluster::spannedCols(const TBCluster* tbcluster, const TBEvent* event)
{
	DUT* dut = event->dut;
	int min = dut->getNcols()-1;
	int max = 0;
	for(auto hit: tbcluster->hits)
	{
		if(hit->col < min)
		{
			min = hit->col;
		}
		if(hit->col > max)
		{
			max = hit->col;
		}
	}
	return max - min + 1;
}

// TODO 
double TBCluster::getEtaCorrectedCol(const TBCluster* tbcluster, const TBEvent* event)
{
	DUT* dut = event->dut;
	double weighted = TBCluster::getChargeWeightedX(tbcluster, event);
	double floored = std::floor(weighted);
	if(floored ==  weighted)
	{
		return weighted;
	}
	return floored + dut->ecorrs.getX(weighted - floored);
}

double TBCluster::getEtaCorrectedX(const TBCluster* tbcluster, const TBEvent* event)
{
	DUT* dut = event->dut;
	double doubleCol = TBCluster::getEtaCorrectedCol(tbcluster, event);
	double doubleRow = TBCluster::getEtaCorrectedRow(tbcluster, event);
	
	double centerX = dut->pixelArray[(int) std::floor(doubleCol)][(int) std::floor(doubleRow)].getLowerLeftX() + doubleCol;
	
	return centerX;
}


double TBCluster::getEtaCorrectedRow(const TBCluster* tbcluster, const TBEvent* event)
{
	DUT* dut = event->dut;
	double weighted = TBCluster::getChargeWeightedY(tbcluster, event);
	double floored = floor(weighted);
	if(floored ==  weighted)
	{
		return weighted;
	}
	return floored + dut->ecorrs.getY(weighted - floored);
}


double TBCluster::getEtaCorrectedY(const TBCluster* tbcluster, const TBEvent* event)
{
	DUT* dut = event->dut;
	double doubleCol = TBCluster::getEtaCorrectedCol(tbcluster, event);
	double doubleRow = TBCluster::getEtaCorrectedRow(tbcluster, event);
	
	double centerY = dut->pixelArray[(int) std::floor(doubleCol)][(int) std::floor(doubleRow)].getLowerLeftY() + doubleRow;
	
	return centerY;
}
// TODO end