#include <fstream>

// ROOT
#include "TMath.h"

// Gaudi
#include "GaudiKernel/PhysicalConstants.h"
#include "GaudiUtils/HistoLabels.h"

// Tb/TbKernel
#include "TbKernel/TbConstants.h"
#include "TbKernel/TbModule.h"

// Local
#include "TbEffPur.h"

using namespace Gaudi::Utils::Histos;

DECLARE_ALGORITHM_FACTORY(TbEffPur)

//=============================================================================
// Standard constructor
//=============================================================================
TbEffPur::TbEffPur(const std::string& name, ISvcLocator* pSvcLocator)
  : TbAlgorithm(name, pSvcLocator),
  m_pitch(0.055),
  m_nDUTpixels(256),
  m_nTracks(0),
  m_nClusters(0),
  m_nTrackedClusters(0),
	m_nClustersPassedCentral(0),
	m_nTracksCentral(0),
	m_nClustersPassedCorner(0),
	m_nTracksCorner(0),
  m_eff(0.),
  m_pur(0.),
  m_deadAreaRadius(2.),
  m_trackAssociated(NULL)
{
  declareProperty("TrackLocation",
                  m_trackLocation = LHCb::TbTrackLocation::Default);
  declareProperty("VertexLocation",
                  m_vertexLocation = LHCb::TbVertexLocation::Default);
  declareProperty("ClusterLocation",
                  m_clusterLocation = LHCb::TbClusterLocation::Default);

  // Parameters.
  declareProperty("DUTindex", m_DUTindex = 4);
  declareProperty("GlobalCutXLow", m_xLow = 3);
  declareProperty("GlobalCutXUp", m_xUp = 9);
  declareProperty("GlobalCutYLow", m_yLow = 3);
  declareProperty("GlobalCutYUp", m_yUp = 9);
  declareProperty("LocalProbabilityCut", m_probCut = 0.5);

  declareProperty("RResidualCut", m_rResidualCut = 0.1);
  declareProperty("TResidualCut", m_tResidualCut = 100);
  declareProperty("ChargeCutLow", m_chargeCutLow = 0);
  // ... note this global cut enforced at different point to above.
  declareProperty("ChargeCutUp", m_chargeCutUp = 1000000);
  declareProperty("LitSquareSide", m_litSquareSide = 0);

  declareProperty("ViewerOutput", m_viewerOutput = false);
  declareProperty("ViewerEventNum", m_viewerEvent = 100);
  declareProperty("TGap", m_tGap = 200);
  declareProperty("CorrelationTimeWindow", m_correlationTimeWindow = 100);
  declareProperty("ApplyVeto", m_applyVeto = true);
  declareProperty("TelescopeClusterVetoDelT", m_telescopeClusterVetoDelT = 30);
  declareProperty("EdgeVetoDistance", m_edgeVetoDistance = 0.025);
  m_nEvent = -1;
}

//=============================================================================
// Initialization
//=============================================================================
StatusCode TbEffPur::initialize() {

  // Initialise the base class.
  StatusCode sc = TbAlgorithm::initialize();
  if (sc.isFailure()) return sc;

	// Efficiency objects.
	m_effX = new TEfficiency("efficiencyX", "efficiencyX", 40, 0.0, 0.11);
	m_effY = new TEfficiency("efficiencyY", "efficiencyY", 40, 0.0, 0.11);
	m_effs = new TEfficiency("efficiency", "efficiency", 1, -0.5, 0.5);
  m_purs = new TEfficiency("purity", "Tb/TbEffPur/pur", 1, -0.5, 0.5);
  m_effHitmap = new TEfficiency("efficiencyHitmap", "efficiencyHitmap",
  		64, 0.0, 14.08, 64, 0.0, 14.08);
  m_purHitmap = new TEfficiency("puritHitmap", "purityHitmap",
  		64, 0.0, 14.08, 64, 0.0, 14.08);

  int nBins = 50;
  m_effHitmapInterPixel = new TEfficiency("efficiencyHitmapInterPixel", "efficiencyHitmapInterPixel",
  		nBins, 0.0, 0.055, nBins, 0.0, 0.055);
  m_purHitmapInterPixel = new TEfficiency("puritHitmapInterPixel", "purityHitmapInterPixel",
  		nBins, 0.0, 0.055, nBins, 0.0, 0.055);

  m_effHitmapInterPixelTriple = new TEfficiency("efficiencyHitmapInterPixelTriple", "efficiencyHitmapInterPixelTriple",
  		150, 0.0, 10*0.055, 20, 0.0, 2*0.055);

  for (unsigned int i=1; i<5; i++) {
  	std::string name = "efficiencyHitmapInterClusterSize" + std::to_string(i);
  	TEfficiency * e = new TEfficiency(name.c_str(), name.c_str(),
  		nBins, 0.0, 0.055, nBins, 0.0, 0.055);
  	m_effHitmapInterPixelVsSizes.push_back(e);
  }

  // Plots for post correlations.
  m_remainsCorrelationsX = book2D("remainsClusterCorrelationsX", "remainsClusterCorrelationsX", 0, 14, 200, 0, 14, 200);
  m_remainsCorrelationsY = book2D("remainsClusterCorrelationsY", "remainsClusterCorrelationsY", 0, 14, 200, 0, 14, 200);
  m_remainsDifferencesXY = book2D("remainsDifferencesXY", "remainsDifferencesXY", -1.5, 1.5, 150, -1.5, 1.5, 150);
  m_clusterRemainsPositionsLocal = book2D("clusterRemainsPositionsLocal", "clusterRemainsPositionsLocal", -5, 20, 600, -5, 20, 600);
  m_trackRemainsPositionsLocal = book2D("trackRemainsPositionsLocal", "trackRemainsPositionsLocal", 0, 256, 600, 0, 256, 600);
  m_clusterRemainsPositionsGlobal = book2D("clusterRemainsPositionsGlobal", "clusterRemainsPositionsGlobal", -5, 20, 600, -5, 20, 600);
  m_trackRemainsPositionsGlobal = book2D("trackRemainsPositionsGlobal", "trackRemainsPositionsGlobal", -5, 20, 600, -5, 20, 600);
  m_vetoTracksHitmap = book2D("vetoTrackHitmap", "vetoTrackHitmap", -5, 20, 600, -5, 20, 600);
  m_vetoClustersHitmap = book2D("vetoClusterHitmap", "vetoClusterHitmap", -5, 20, 600, -5, 20, 600);
  m_timeResidualVsColumn = book2D("timeResidualVsColumn", "timeResidualVsColumn", 0, 256, 256, -50, 50, 500);


  // Setup the tools.
  m_trackFit = tool<ITbTrackFit>("TbTrackFit", "Fitter", this);
  // h_effVsCharge = new TH2F("effVsCharge", "effVsCharge", 50, 0.5, 1, 20, 0.5, 20.5);
  // TFile * f = new TFile("/afs/cern.ch/user/d/dsaunder/cmtuser/KEPLER/KEPLER_HEAD/EDR_plots/EfficiencyResults/Eff-S22-Run6309-500V.root", "READ");
  // h_effInter = (TH2F*)f->Get("eff_histo");
  return StatusCode::SUCCESS;

}


//=============================================================================
// Finalizer
//=============================================================================
StatusCode TbEffPur::finalize() {
	m_effs->SetTotalEvents(0, m_nTracks);
	m_effs->SetPassedEvents(0, m_nTrackedClusters);

	m_purs->SetTotalEvents(0, m_nClusters);
	m_purs->SetPassedEvents(0, m_nTrackedClusters);

	m_eff = m_nTrackedClusters/(1.0*m_nTracks);
	m_pur = m_nTrackedClusters/(1.0*m_nClusters);
	std::cout<<"ith ROOT Efficiency: " << m_effs->GetEfficiency(0) << " + " << m_effs->GetEfficiencyErrorUp(0) <<" - "<< m_effs->GetEfficiencyErrorLow(0) <<std::endl;
	std::cout<<"ith ROOT Purity: " << m_purs->GetEfficiency(0) << " + " << m_purs->GetEfficiencyErrorUp(0) <<" - "<< m_purs->GetEfficiencyErrorLow(0) <<std::endl;

	std::cout<<"Corner efficiency: "<<m_nClustersPassedCorner/(1.*m_nTracksCorner)<<std::endl;
	std::cout<<"Central efficiency: "<<m_nClustersPassedCentral/(1.*m_nTracksCentral)<<std::endl;

	plot(m_eff, "Efficiency", "Efficiency", 0.0, 1, 1);
	plot(m_pur, "Purity", "Purity", 0.0, 1, 1);
	plot(m_effs->GetEfficiencyErrorUp(0), "EfficiencyError", "EfficiencyError", 0.0, 1, 1);
	plot(m_purs->GetEfficiencyErrorUp(0), "PurityError", "PurityError", 0.0, 1, 1);

	TH2D * h = Gaudi::Utils::Aida2ROOT::aida2root(m_trackRemainsPositionsLocal);
	double maxBin = h->GetBinContent(h->GetMaximumBin());
	double minBin = h->GetBinContent(h->GetMinimumBin());
	for (int i=0; i<m_trackRemainsPositionsLocal->xAxis().bins(); i++) {
		for (int j=0; j<m_trackRemainsPositionsLocal->yAxis().bins(); j++) {
			plot(m_trackRemainsPositionsLocal->binHeight(i, j), "trackRemainsBins",
					"trackRemainsBins", minBin, maxBin, int(maxBin-minBin));
		}
	}

	TFile * f = new TFile("EfficiencyResults.root", "RECREATE");
	//m_effs->Write();

	m_effHitmapInterPixel->Write();
	m_effHitmapInterPixelTriple->Write();
	//h_effVsCharge->Write();
	//m_effs->CreateHistogram()->Write();
	m_effHitmapInterPixel->CreateHistogram()->Write();
	m_effHitmapInterPixelTriple->CreateHistogram()->Write();



//	for (unsigned int i=0; i<m_effHitmapInterPixelVsSizes.size(); i++) {
//		m_effHitmapInterPixelVsSizes[i]->Write();
//		//m_effHitmapInterPixelVsSizes[i]->CreateHistogram()->Write();
//	}
	f->Close();
  return TbAlgorithm::finalize();
}


//=============================================================================
// Main execution
//=============================================================================
StatusCode TbEffPur::execute()
{
	m_nEvent++;

  m_tracks = getIfExists<LHCb::TbTracks>(m_trackLocation);
  //m_vertices = getIfExists<LHCb::TbVertices>(m_vertexLocation);
  m_clusters = getIfExists<LHCb::TbClusters>(m_clusterLocation + std::to_string(m_DUTindex));

  if (m_clusters->size() < 100 || m_tracks->size() < 100) return StatusCode::SUCCESS;

  // Matching.
  effPur();
  return StatusCode::SUCCESS;
}


//=============================================================================
// Efficiency finding
//=============================================================================
void TbEffPur::effPur()
{
	// Non-veto'd containers.
  std::vector<LHCb::TbCluster*> * cutClusters = new std::vector<LHCb::TbCluster*>();
  std::vector<LHCb::TbTrack*> cutTracks;
  int nClustersPreVeto = m_clusters->size();
  int nTracksPreVeto  = m_tracks->size();


  // Apply veto.
  if (m_applyVeto) applyVeto(cutClusters, &cutTracks);
  else fillAllTrackClusters(cutClusters, &cutTracks);
  m_trackAssociated = new std::vector<bool>(cutTracks.size(), false);
  //std::cout<<"Fraction tracks not veto'd: "<<cutTracks.size()/(1.0*nTracksPreVeto)<<"\tEvent: "<<m_nEvent<<std::endl;
  //std::cout<<"Fraction clusters not veto'd: "<<cutClusters->size()/(1.0*nClustersPreVeto)<<"\tEvent: "<<m_nEvent<<std::endl;
  plot(cutClusters->size()/(1.0*nClustersPreVeto), "clusterFractionVetod", "clusterFractionVetod", 0.0, 1.0, 200);
  plot(cutTracks.size()/(1.0*nTracksPreVeto), "trackFractionVetod", "trackFractionVetod", 0.0, 1.0, 200);


  // Do matching.
  m_nClusters += cutClusters->size();
  m_nTracks += cutTracks.size();
  trackClusters(cutClusters, &cutTracks);


  // Viewer options.
  if (m_nEvent == m_viewerEvent) {
  	if (m_viewerOutput) outputViewerData();
		for (LHCb::TbClusters::iterator iclust = cutClusters->begin();
				iclust != cutClusters->end(); iclust++) {
//			if (!(*iclust)->associated())
//				std::cout<<"Note: non-associated clusters at hTime: "<<(*iclust)->htime()<<std::endl;
		}
  }

  // Correlate remains.
  //correlateRemains(&cutTracks, cutClusters);
  delete cutClusters;
  delete m_trackAssociated;
}


//=============================================================================
// Global cuts
//=============================================================================
void TbEffPur::applyVeto(std::vector<LHCb::TbCluster*> * cutClusters,
		std::vector<LHCb::TbTrack*> * cutTracks)
{
  // Gather all cluster times (over all planes) and sort it ___________________
	std::vector<double> allClusterTimes;
	for (unsigned int i=0; i<m_nPlanes; i++) {
    LHCb::TbClusters * clusters = getIfExists<LHCb::TbClusters>(m_clusterLocation + std::to_string(i));
		for (LHCb::TbClusters::const_iterator it = clusters->begin(); it != clusters->end(); ++it)
			allClusterTimes.push_back((*it)->htime());
	}
	std::sort(allClusterTimes.begin(), allClusterTimes.end());



	// Get big enough gaps in pairs _____________________________________________
	std::vector<double> tGapCenters;
	tGapCenters.push_back(0);
	for (std::vector<double>::iterator itime = allClusterTimes.begin()+1;
			itime != allClusterTimes.end(); itime++) {
		plot((*itime) - (*(itime-1)), "allGapWidths", "allGapWidths", 0.0, 3000.0, 200);
		if ((*itime) - (*(itime-1)) > m_tGap) {
			tGapCenters.push_back(0.5*((*itime) + (*(itime-1))));
			plot((*itime) - (*(itime-1)), "bigGapWidths", "bigGapWidths", 0.0, 3000.0, 200);
		}
	}
	tGapCenters.push_back(allClusterTimes.back());
	counter("nGaps") += tGapCenters.size();
	plot(tGapCenters.size(), "nGaps", "nGaps", 0.0, 1000.0, 200);



	// Loop over these sub events and veto ______________________________________
	std::vector<double>::iterator igap;
	for (igap = tGapCenters.begin(); igap != tGapCenters.end() - 1; igap++) {
		bool veto = false;

		if (igap == tGapCenters.begin() || igap == tGapCenters.end() -1) veto = true;

		// Loop tracks.
		if (!veto) {
			for (LHCb::TbTrack* track : *m_tracks) {
				if (track->htime() > (*igap) && track->htime() < (*(igap+1))) {
					// Veto on tracks outside cuts.
					const Gaudi::XYZPoint trackInterceptGlobal = geomSvc()->intercept(track, m_DUTindex);
					if (!globalCutPosition(trackInterceptGlobal)) veto = true;
					else if (outsideDUT(trackInterceptGlobal)) veto = true;
					else if (interceptDeadPixel(trackInterceptGlobal)) veto = true;
				}
			}
		}


		// Loop vertices.
//		if (!veto) {
//			for (LHCb::TbVertex* vertex : *m_vertices) {
//				if (vertex->htime() > (*igap) && vertex->htime() < (*(igap+1))) {
//					veto = true;
//				}
//			}
//		}


		// Loop DUT clusters.
		if (!veto) {
			for (const LHCb::TbCluster* cluster : *m_clusters) {
				if (cluster->htime() > (*igap) && cluster->htime() < (*(igap+1))) {
					// Veto on clusters outside cuts.
					Gaudi::XYZPoint clusterInterceptGlobal(cluster->x(), cluster->y(), cluster->z());
					if (!globalCutPosition(clusterInterceptGlobal)) veto = true;
					if (cluster->charge() > m_chargeCutUp || cluster->charge() < m_chargeCutLow) veto = true;
				}
			}
		}


		// Loop telescope clusters.
		if (!veto) {
			for (unsigned int i=0; i<m_nPlanes; i++) {
				if (i == m_DUTindex) continue;
				LHCb::TbClusters * clusters = getIfExists<LHCb::TbClusters>(m_clusterLocation + std::to_string(i));
				for (LHCb::TbClusters::const_iterator it = clusters->begin(); it != clusters->end()-1; ++it) {
					if ((*it)->htime() > (*igap) && (*it)->htime() < (*(igap+1))) {
						double delt = (*it)->htime() - (*(it+1))->htime();
						if (fabs(delt) < m_telescopeClusterVetoDelT) {
							veto = true;
							break;
						}
					}
				}
			}
		}


		if (!veto) fillTrackClusters(cutClusters, cutTracks, (*igap), (*(igap+1)));
		else counter("nGapsFiltered") += 1;
	}
}



//=============================================================================
// Correlate remains
//=============================================================================
void TbEffPur::correlateRemains(std::vector<LHCb::TbTrack*> * cutTracks,
		std::vector<LHCb::TbCluster*> * cutClusters)
{
	for (std::vector<LHCb::TbCluster*>::iterator ic = cutClusters->begin();
	  		ic != cutClusters->end(); ++ic) {
		if (!(*ic)->associated()) {
			m_clusterRemainsPositionsGlobal->fill((*ic)->x(), (*ic)->y());
			m_clusterRemainsPositionsLocal->fill((*ic)->xloc(), (*ic)->yloc());
			plot((*ic)->charge(), "chargeClusterRemains", "chargeClusterRemains", 0.0, 1000.0, 125);
		}
	}

	unsigned int i = 0;
	for (std::vector<LHCb::TbTrack*>::iterator itrack = cutTracks->begin();
			itrack != cutTracks->end(); itrack++) {
		if (!m_trackAssociated->at(i)) {
		  Gaudi::XYZPoint trackIntercept = geomSvc()->intercept((*itrack), m_DUTindex);
		  m_trackRemainsPositionsGlobal->fill(trackIntercept.x(), trackIntercept.y());
		  auto interceptUL = geomSvc()->globalToLocal(trackIntercept, m_DUTindex);
		  m_trackRemainsPositionsLocal->fill(interceptUL.x()/m_pitch - 0.5, interceptUL.y()/m_pitch - 0.5);
		}
		i++;
	}


  // Draws correlations plots between non-associated clusters and tracks.
	for (std::vector<LHCb::TbTrack*>::iterator itrack = cutTracks->begin();
			itrack != cutTracks->end(); itrack++) {
	  for (std::vector<LHCb::TbCluster*>::iterator ic = cutClusters->begin();
	  		ic != cutClusters->end(); ++ic) {
	  	if ((*ic)->htime() < (*itrack)->htime() - m_correlationTimeWindow) continue;
	  	if (!(*ic)->associated()) {
	  		//if (m_nEvent == m_viewerEvent) std::cout<<"Impurity at time: "<<(*ic)->htime()<<std::endl;
	  		Gaudi::XYZPoint trackIntercept = geomSvc()->intercept((*itrack), m_DUTindex);
	  		m_remainsCorrelationsX->fill(trackIntercept.x(), (*ic)->x());
	  		m_remainsCorrelationsY->fill(trackIntercept.y(), (*ic)->y());

	  		plot(trackIntercept.x() - (*ic)->x(), "remainClustDifferencesX", "remainsClustDifferencesX", -1.5, 1.5, 150);
	  		plot(trackIntercept.y() - (*ic)->y(), "remainClustDifferencesY", "remainClustDifferencesY", -1.5, 1.5, 150);
	  		plot((*itrack)->htime() - (*ic)->htime(), "remainClustDifferencesT", "remainClustDifferencesT", -50, 50, 150);

	  		m_remainsDifferencesXY->fill(trackIntercept.y() - (*ic)->y(), trackIntercept.x() - (*ic)->x());
	  	}
	  	if ((*ic)->htime() > (*itrack)->htime() + m_correlationTimeWindow) break;
	  }
	}
}


//=============================================================================
//
//=============================================================================
void TbEffPur::fillAllTrackClusters(std::vector<LHCb::TbCluster*> * cutClusters,
		std::vector<LHCb::TbTrack*> * cutTracks)
{
  for (LHCb::TbTrack* track : *m_tracks) cutTracks->push_back(track);
  for (std::vector<LHCb::TbCluster*>::iterator iClust = m_clusters->begin();
  		iClust != m_clusters->end(); iClust++) cutClusters->push_back(*iClust);
}


//=============================================================================
//
//=============================================================================
void TbEffPur::fillTrackClusters(std::vector<LHCb::TbCluster*> * cutClusters,
		std::vector<LHCb::TbTrack*> * cutTracks, double tlow, double tup)
{
	// Push tracks and clusters inside this time window, which passed the veto.
  for (LHCb::TbTrack* track : *m_tracks) {
  	if (track->htime() > tlow && track->htime() < tup) {
  		Gaudi::XYZPoint trackIntercept = geomSvc()->intercept(track, m_DUTindex);
  		cutTracks->push_back(track);
  		m_vetoTracksHitmap->fill(trackIntercept.x(), trackIntercept.y());

  		Gaudi::XYZPoint trackInterceptLocal = geomSvc()->globalToLocal(trackIntercept, m_DUTindex);
			int row = trackInterceptLocal.y()/m_pitch - 0.5;
			int col = trackInterceptLocal.x()/m_pitch - 0.5;

			if (pixelSvc()->isMasked(pixelSvc()->address(col, row), m_DUTindex)) {
				std::cout<<"Shouldn't be here!"<<std::endl;
			}
  	}
  }

  for (LHCb::TbCluster* cluster : *m_clusters) {
		if (cluster->htime() > tlow && cluster->htime() < tup) {
			cutClusters->push_back(cluster);
			m_vetoClustersHitmap->fill(cluster->x(), cluster->y());
		}
  }
}


//=============================================================================
// Global cuts
//=============================================================================
void TbEffPur::trackClusters(std::vector<LHCb::TbCluster*> * cutClusters,
		std::vector<LHCb::TbTrack*> * cutTracks)
{
  for (std::vector<LHCb::TbTrack*>::iterator itrack = cutTracks->begin();
  		itrack != cutTracks->end(); itrack++) {

  	// Get local position of track.
  	const auto interceptUG = geomSvc()->intercept((*itrack), m_DUTindex);
    const auto interceptUL = geomSvc()->globalToLocal(interceptUG, m_DUTindex);


  	// Loop over clusters to find the closest.
  	LHCb::TbCluster * closestCluster = NULL;
  	double bestRadialDistance;
  	for (std::vector<LHCb::TbCluster*>::iterator iclust = cutClusters->begin();
  			iclust != cutClusters->end(); iclust++) {
  		bool match = matchTrackToCluster((*iclust), (*itrack));
  		double radialDistance = getRadialSeparation((*iclust), (*itrack));
  		if (match) {
  			if (!closestCluster) {
  				closestCluster = (*iclust);
  				bestRadialDistance = radialDistance;
  			}
  			else if (radialDistance<bestRadialDistance) {
  				closestCluster = (*iclust);
  				bestRadialDistance = radialDistance;
  			}
  		}

  		if ((*iclust)->htime() - (*itrack)->htime() > m_tResidualCut) break;
  	}

  	double trackX = interceptUL.x();
  	double trackY = interceptUL.y();

		if (closestCluster != NULL && closestCluster->charge() > m_chargeCutLow) {
			closestCluster->setAssociated(true);
			m_nTrackedClusters++;
			m_effHitmap->Fill(true, trackX, trackY);
			m_effX->Fill(true, trackX);
			m_effY->Fill(true, trackY);

			if (fmod(trackX, m_pitch)<0.0183 && fmod(trackY, m_pitch) < 0.0183) m_nClustersPassedCorner++;
			if (fmod(trackX, m_pitch) > 0.0183 &&
					fmod(trackX, m_pitch) < (2*0.0183) &&
					fmod(trackY, m_pitch) > 0.0183 &&
					fmod(trackY, m_pitch) < (2*0.0183)) m_nClustersPassedCentral++;


			//h_effVsCharge->Fill(h_effInter->Interpolate(fmod(trackX, m_pitch), fmod(trackY, m_pitch)), closestCluster->charge());
			// if (trackY > 9.25 && trackY < 10.25) m_effHitmapInterPixel->Fill(true, fmod(trackX, m_pitch), fmod(trackY, m_pitch));
			// if (trackX > 28.05 &&
			// 		trackX < 28.6) m_effHitmapInterPixelTriple->Fill(true, trackX-28.05, fmod(trackY, m_pitch));


      m_effHitmapInterPixel->Fill(true, fmod(trackX, m_pitch), fmod(trackY, m_pitch));
      m_effHitmapInterPixelTriple->Fill(true, trackX-28.05, fmod(trackY, m_pitch));


			if (closestCluster->size() <= m_effHitmapInterPixelVsSizes.size())
				m_effHitmapInterPixelVsSizes[closestCluster->size()-1]->Fill(true, fmod(trackX, m_pitch), fmod(trackY, m_pitch));
			m_trackAssociated->at(itrack - cutTracks->begin()) = true;
			plot(closestCluster->x() - interceptUG.x(), "xResidualsClustMinusTrack", "xResidualsClustMinusTrack", -0.2, 0.2, 400);
			plot(closestCluster->y() - interceptUG.y(), "yResidualsClustMinusTrack", "yResidualsClustMinusTrack", -0.2, 0.2, 400);
			plot(closestCluster->htime() - (*itrack)->htime(), "hTimeResidualClustMinusTrack", "hTimeResidualClustMinusTrack", -50, 50, 400);
			if (closestCluster->size() == 1) {
				auto hits = closestCluster->hits();
				int col = (*hits.begin())->col();
				m_timeResidualVsColumn->fill(col, closestCluster->htime() - (*itrack)->htime());
			}
			if ((*itrack)->vertexed()) counter("nVerticesCorrelated") += 1;
		}
		else {
			m_effHitmap->Fill(false, trackX, trackY);
			// if (trackY > 9.25 && trackY < 10.25) m_effHitmapInterPixel->Fill(false, fmod(trackX, m_pitch), fmod(trackY, m_pitch));
			// if (trackX > 28.05 &&
			// 		trackX < 28.6) m_effHitmapInterPixelTriple->Fill(false, trackX-28.05, fmod(trackY, m_pitch));

      m_effHitmapInterPixel->Fill(false, fmod(trackX, m_pitch), fmod(trackY, m_pitch));
      m_effHitmapInterPixelTriple->Fill(false, trackX-28.05, fmod(trackY, m_pitch));
			m_effX->Fill(false, trackX);
			m_effY->Fill(false, trackY);
			//std::cout<<"Inefficiency at time: "<<(*itrack)->htime()<<std::endl;
		}
		if (fmod(trackX, m_pitch)<0.0183 && fmod(trackY, m_pitch) < 0.0183) m_nTracksCorner++;
		if (fmod(trackX, m_pitch) > 0.0183 &&
					fmod(trackX, m_pitch) < (2*0.0183) &&
					fmod(trackY, m_pitch) > 0.0183 &&
					fmod(trackY, m_pitch) < (2*0.0183)) m_nTracksCentral++;
  }
}



//=============================================================================
//
//=============================================================================
double TbEffPur::getRadialSeparation(LHCb::TbCluster * cluster,
		LHCb::TbTrack * track)
{
  Gaudi::XYZPoint trackIntercept = geomSvc()->intercept(track, m_DUTindex);
	double xResidual = cluster->x() - trackIntercept.x();
	double yResidual = cluster->y() - trackIntercept.y();

	double r2 = xResidual*xResidual + yResidual*yResidual;
	return pow(r2, 0.5);
}


//=============================================================================
// Local cuts
//=============================================================================
bool TbEffPur::matchTrackToCluster(LHCb::TbCluster * cluster,
		LHCb::TbTrack * track)
{
  Gaudi::XYZPoint trackIntercept = geomSvc()->intercept(track, m_DUTindex);
	double xResidual = cluster->x() - trackIntercept.x();
	double yResidual = cluster->y() - trackIntercept.y();
	double tResidual = cluster->htime() - track->htime();
	double delr = pow(pow(xResidual, 2) + pow(yResidual, 2), 0.5);

//	if (track->vertexed()) plot(xResidual, "vertexResiduals/X", "vertexResiduals/X", -0.5, 0.5, 200);
//	if (track->vertexed()) plot(yResidual, "vertexResiduals/Y", "vertexResiduals/Y", -0.5, 0.5, 200);
//	if (track->vertexed()) plot(tResidual, "vertexResiduals/T", "vertexResiduals/T", -0.5, 0.5, 200);
//
//	if (track->vertexed() && litPixel(cluster, track)) {
//		plot(xResidual, "vertexResiduals/Xlit", "vertexResiduals/Xlit", -0.5, 0.5, 200);
//		plot(yResidual, "vertexResiduals/Ylit", "vertexResiduals/Ylit", -0.5, 0.5, 200);
//	}


	if (fabs(tResidual) > m_tResidualCut) {
  	counter("nTimeRejected") += 1;
		return false;
	}

  if (delr > m_rResidualCut) {
  	counter("nSpatialRejected") += 1;
  	if (litPixel(cluster, track)) return true;
  	else {
  		counter("nSpatialAndLitRejected") += 1;
  		return false;
  	}
  }

//	if (!litPixel(cluster, track)) return false;

  return true;
}


//=============================================================================
//
//=============================================================================
bool TbEffPur::interceptDeadPixel(Gaudi::XYZPoint trackInterceptGlobal)
{
	// Find the row and column corresponding to the track intercept.
	Gaudi::XYZPoint trackInterceptLocal = geomSvc()->globalToLocal(trackInterceptGlobal, m_DUTindex);
  int row = trackInterceptLocal.y()/m_pitch - 0.5;
  int col = trackInterceptLocal.x()/m_pitch - 0.5;

  for (int icol = col - 1; icol != col+2; icol++) {
  	for (int irow = row - 1; irow != row+2; irow++) {
  		if (icol >= 0 && icol <256 && irow>=0 && irow<256) {
  			if (pixelSvc()->isMasked(pixelSvc()->address(icol, irow), m_DUTindex)) {
					return true;
				}

//  			if (icol == 17 && irow == 36) {
//  				return true;
//  			}
//				if (icol == 18 && irow == 36) {
//  				return true;
//  			}
//				if (icol == 53 && irow == 3) {
//  				return true;
//  			}
//				if (icol == 54 && irow == 3) {
//  				return true;
//  			}
//				if (icol == 73 && irow == 26) {
//  				return true;
//  			}
//				if (icol == 74 && irow == 26) {
//  				return true;
//  			}
//				if (icol == 19 && irow == 103) {
//  				return true;
//  			}
//				if (icol == 31 && irow == 106) {
//  				return true;
//  			}
//				if (icol == 32 && irow == 106) {
//  				return true;
//  			}
//				if (icol == 38 && irow == 108) {
//  				return true;
//  			}
//				if (icol == 63 && irow == 95) {
//  				return true;
//  			}
//				if (icol == 64 && irow == 95) {
//  				return true;
//  			}
//				if (icol == 103 && irow == 23) {
//  				return true;
//  			}
//				if (icol == 104 && irow == 23) {
//  				return true;
//  			}
//				if (icol == 115 && irow == 20) {
//  				return true;
//  			}
//				if (icol == 116 && irow == 94) {
//  				return true;
//  			}
  		}
  	}
  }



  return false;
}

//=============================================================================
//
//=============================================================================
bool TbEffPur::litPixel(LHCb::TbCluster * cluster, LHCb::TbTrack * track)
{
	// Find the row and column corresponding to the track intercept.
	Gaudi::XYZPoint trackIntercept = geomSvc()->intercept(track, m_DUTindex);
	Gaudi::XYZPoint pLocal = geomSvc()->globalToLocal(trackIntercept, m_DUTindex);
  int row = pLocal.y()/m_pitch - 0.5;
  int col = pLocal.x()/m_pitch - 0.5;

  // See if within a pre-defined square of pixels.
  for (unsigned int i=0; i<cluster->size(); i++) {
  	unsigned int delRow = abs(row - int(cluster->hits()[i]->row()));
  	unsigned int delCol = abs(col - int(cluster->hits()[i]->col()));

  	if (delRow <= m_litSquareSide && delCol <= m_litSquareSide) {
		  counter("nLitPixels") += 1;
  		return true;
  	}
  }

  return false;
}


//=============================================================================
// Excluding dead regions.
//=============================================================================
bool TbEffPur::globalCutPosition(Gaudi::XYZPoint pGlobal)
{
	if (pGlobal.x() < m_xLow + m_edgeVetoDistance) return false;
	else if (pGlobal.x() > m_xUp - m_edgeVetoDistance) return false;
	else if (pGlobal.y() < m_yLow + m_edgeVetoDistance) return false;
	else if	(pGlobal.y() > m_yUp - m_edgeVetoDistance) return false;
  return true; // Inside.
}


//=============================================================================
// Excluding dead regions.
//=============================================================================
bool TbEffPur::outsideDUT(Gaudi::XYZPoint interceptUG)
{
	const auto interceptUL = geomSvc()->globalToLocal(interceptUG, m_DUTindex);
	if (interceptUL.x() < 0 + m_edgeVetoDistance ||
			interceptUL.x() > 14.08 - m_edgeVetoDistance ||
			interceptUL.y() < 0 + m_edgeVetoDistance||
			interceptUL.y() > 14.08 - m_edgeVetoDistance) return true;
	return false;
}



//=============================================================================
// Viewer output.
//=============================================================================
//=============================================================================
// Viewer output.
//=============================================================================
//=============================================================================
// Viewer output.
//=============================================================================

void TbEffPur::outputDeadRegion(unsigned int col, unsigned int row) {
	std::ofstream myfile;
  myfile.open("/afs/cern.ch/user/d/dsaunder/KeplerViewerData.dat", std::ofstream::app);
	std::string outputLine;
	for (int irow = std::max(int(row-m_deadAreaRadius), 0);
  		irow < std::min(int(row+m_deadAreaRadius+1), m_nDUTpixels); irow++) {
  	for (int icol = std::max(int(col-m_deadAreaRadius), 0);
  		icol < std::min(int(col+m_deadAreaRadius+1), m_nDUTpixels); icol++) {


			outputLine = "DeadPixel ";
			const double xLocal = (col-icol) * m_pitch; // Dont want the middle!
			const double yLocal = (row-irow) * m_pitch; // Dont want the middle!
			Gaudi::XYZPoint pLocal(xLocal, yLocal, 0.);

			Gaudi::XYZPoint posn = geomSvc()->localToGlobal(pLocal, m_DUTindex);
			outputLine += std::to_string(posn.x()) + " ";
			outputLine += std::to_string(posn.y()) + " ";
			outputLine += std::to_string(posn.z()) + " ";

			Gaudi::XYZPoint posn2(pLocal.x() + 0.055, pLocal.y(), 0.);
			posn = geomSvc()->localToGlobal(posn2, m_DUTindex);
			outputLine += std::to_string(posn.x()) + " ";
			outputLine += std::to_string(posn.y()) + " ";
			outputLine += std::to_string(posn.z()) + " ";

			Gaudi::XYZPoint posn3(pLocal.x() + 0.055, pLocal.y()+0.055, 0.);
			posn = geomSvc()->localToGlobal(posn3, m_DUTindex);
			outputLine += std::to_string(posn.x()) + " ";
			outputLine += std::to_string(posn.y()) + " ";
			outputLine += std::to_string(posn.z()) + " ";

			Gaudi::XYZPoint posn4(pLocal.x(), pLocal.y()+0.055, 0.);
			posn = geomSvc()->localToGlobal(posn4, m_DUTindex);
			outputLine += std::to_string(posn.x()) + " ";
			outputLine += std::to_string(posn.y()) + " ";
			outputLine += std::to_string(posn.z()) + " ";

			outputLine += "\n";
			myfile << outputLine;
  	}
	}
	myfile.close();
}


//=============================================================================
// Viewer output.
//=============================================================================

void TbEffPur::outputViewerData()
{
  std::ofstream myfile;
  myfile.open("/afs/cern.ch/user/d/dsaunder/KeplerViewerData.dat", std::ofstream::app);
  std::string outputLine;
//
//  // First output the chips.
//  for (unsigned int i=0; i<m_nPlanes; i++) {
//  	outputLine = "Chip ";
//  	Gaudi::XYZPoint posn1(0., 14.08, 0.);
//    Gaudi::XYZPoint posn = geomSvc()->localToGlobal(posn1, i);
//    outputLine += std::to_string(posn.x()) + " ";
//    outputLine += std::to_string(posn.y()) + " ";
//    outputLine += std::to_string(posn.z()) + " ";
//
//    Gaudi::XYZPoint posn2(14.08, 14.08, 0.);
//    posn = geomSvc()->localToGlobal(posn2, i);
//    outputLine += std::to_string(posn.x()) + " ";
//    outputLine += std::to_string(posn.y()) + " ";
//    outputLine += std::to_string(posn.z()) + " ";
//
//    Gaudi::XYZPoint posn3(14.08, 0., 0.);
//    posn = geomSvc()->localToGlobal(posn3, i);
//    outputLine += std::to_string(posn.x()) + " ";
//    outputLine += std::to_string(posn.y()) + " ";
//    outputLine += std::to_string(posn.z()) + " ";
//
//    Gaudi::XYZPoint posn4(0., 0., 0.);
//    posn = geomSvc()->localToGlobal(posn4, i);
//    outputLine += std::to_string(posn.x()) + " ";
//    outputLine += std::to_string(posn.y()) + " ";
//    outputLine += std::to_string(posn.z()) + " ";
//
//    outputLine += "\n";
//    myfile << outputLine;
//  }


  // Clusters.
  auto ic = m_clusters->begin();
	const auto end = m_clusters->end();
	for (; ic != end; ++ic) {
		outputLine = "Cluster ";
		outputLine += std::to_string((*ic)->x()) + " ";
		outputLine += std::to_string((*ic)->y()) + " ";
		outputLine += std::to_string((*ic)->z()) + " ";
		outputLine += std::to_string((*ic)->htime()) + " ";

		// if ((*ic)->endCluster() && (*ic)->vertexed()) outputLine += "5 \n";
		//if ((*ic)->vertexed()) outputLine += "4 \n";
//			else if ((*ic)->endCluster()) outputLine += "3 \n";
		if ((*ic)->associated()) outputLine += "2 \n";
//			else if ((*ic)->volumed()) outputLine += "1 \n";

		else {
			outputLine += "0 \n";
		}
		myfile << outputLine;
	}

	outputLine = "CentralRegion ";
	outputLine += std::to_string(m_xLow) + " ";
	outputLine += std::to_string(m_xUp) + " ";
	outputLine += std::to_string(m_yLow) + " ";
	outputLine += std::to_string(m_yUp) + " ";
	outputLine += std::to_string(geomSvc()->module(m_DUTindex)->z()) + " ";
	myfile << outputLine;

//		// Its hits.
//		for (auto hit : (*ic)->hits()) {
//			outputLine = "Pixel ";
//			const double xLocal = (hit->col()) * m_pitch; // Dont want the middle!
//			const double yLocal = (hit->row()) * m_pitch; // Dont want the middle!
//			Gaudi::XYZPoint pLocal(xLocal, yLocal, 0.);
//
//			Gaudi::XYZPoint posn = geomSvc()->localToGlobal(pLocal, i);
//			outputLine += std::to_string(posn.x()) + " ";
//			outputLine += std::to_string(posn.y()) + " ";
//			outputLine += std::to_string(posn.z()) + " ";
//
//			Gaudi::XYZPoint posn2(pLocal.x() + 0.055, pLocal.y(), 0.);
//			posn = geomSvc()->localToGlobal(posn2, i);
//			outputLine += std::to_string(posn.x()) + " ";
//			outputLine += std::to_string(posn.y()) + " ";
//			outputLine += std::to_string(posn.z()) + " ";
//
//			Gaudi::XYZPoint posn3(pLocal.x() + 0.055, pLocal.y()+0.055, 0.);
//			posn = geomSvc()->localToGlobal(posn3, i);
//			outputLine += std::to_string(posn.x()) + " ";
//			outputLine += std::to_string(posn.y()) + " ";
//			outputLine += std::to_string(posn.z()) + " ";
//
//			Gaudi::XYZPoint posn4(pLocal.x(), pLocal.y()+0.055, 0.);
//			posn = geomSvc()->localToGlobal(posn4, i);
//			outputLine += std::to_string(posn.x()) + " ";
//			outputLine += std::to_string(posn.y()) + " ";
//			outputLine += std::to_string(posn.z()) + " ";
//
//			outputLine += std::to_string(hit->htime()) + " ";
//			outputLine += std::to_string(hit->ToT()) + " ";
//
//			outputLine += "\n";
//			myfile << outputLine;
//		}

  myfile.close();
}


//=============================================================================
// END
//=============================================================================