#include "hstart2clusterfinder.h"

#include <iostream>

#include <assert.h>
#include <algorithm>
#include <fstream>

//_HADES_CLASS_DESCRIPTION
//###############################################################
//###############################################################
// Helper classes for cluster finding in START LGAD
//###############################################################
//###############################################################


//###############################################################
start_singlestriphit::start_singlestriphit()
{
        fmoduleid  = -1;
        fchannelid = -1;
        fmultiplid = -1;
        ftoa       = -1;
        ftot       = -1;
        fstartEntry= -1;
}

start_singlestriphit::start_singlestriphit(Int_t modid,Int_t channel, Int_t multiplid,Double_t toa, Double_t tot,Int_t startEntry)
{
    fmoduleid  = modid;
    fchannelid = channel;
    fmultiplid = multiplid;
    ftoa       = toa;
    ftot       = tot;
    fstartEntry=startEntry;
};
//###############################################################

//###############################################################
start_cluster::start_cluster(std::vector<start_singlestriphit> &hits)
{
    fclusterstrips = hits;
    calccluster();
}

void start_cluster::setclusterstrips(std::vector<start_singlestriphit> &hits)
{
    fclusterstrips = hits;
    calccluster();
}

void start_cluster::calccluster()
{
    assert(fclusterstrips.size() > 0 && "empty cluster");

    Int_t currstripnr = fclusterstrips[0].getModID() * 100 + fclusterstrips[0].getChannelID();
    Int_t maxindex = 0,maxchannel = currstripnr;

    std::vector<Double_t> toarange,totrange;
    std::vector<Int_t>    channelrange;
    for(Int_t i = 0; i < 2; i ++) {
	toarange    .push_back(fclusterstrips[0].getTOA());
	totrange    .push_back(fclusterstrips[0].getTOT());
	channelrange.push_back(currstripnr);
    }

    for(UInt_t i = 0; i < fclusterstrips.size(); i ++) {
	Int_t currstripnr = fclusterstrips[i].getModID() * 100 + fclusterstrips[i].getChannelID();
	// channel range
	if(currstripnr>channelrange[1]) { channelrange[1] = currstripnr; }
	if(currstripnr<channelrange[0]) { channelrange[0] = currstripnr; }

	//calc tot range and max
	if(fclusterstrips[i].getTOT() > totrange[1]) {
	    maxindex    = i;
	    totrange[1] = fclusterstrips[i].getTOT();
	    maxchannel  = currstripnr;
	}
	if(fclusterstrips[i].getTOT() < totrange[0]) { totrange[0] = fclusterstrips[i].getTOT(); }

	// toa range
	if(fclusterstrips[i].getTOA()>toarange[1]) { toarange[1]=fclusterstrips[i].getTOA(); }
	if(fclusterstrips[i].getTOA()<toarange[0]) { toarange[0]=fclusterstrips[i].getTOA(); }
    }
    fToArange       = abs(toarange[1] - toarange[0]);
    fToTrange       = abs(totrange[1] - totrange[0]);
    fchannelrange   = abs(channelrange[1] - channelrange[0]) + 1;
    fmaxsignalstrip = fclusterstrips[maxindex];
    fMaxTot         = totrange[1];
    fMaxstripposincluster = abs(maxchannel - channelrange[0]) + 1; // starts at 1
    fClustersize          = fclusterstrips.size();
}
//###############################################################


//###############################################################
HStart2ClusterFinder::HStart2ClusterFinder()
{
    fToAthresh      = 2;
    fStriprange     = 5;
    fMaxClustersize = 4;
    fclustermaxstrips.resize(0);
    fclustersizes    .resize(0);
}

void HStart2ClusterFinder::setclusterthresholds(Double_t toarange,Int_t striprange, UInt_t maxclustersize)
{
    fToAthresh      = toarange;
    fStriprange     = striprange;
    fMaxClustersize = maxclustersize;
}

std::vector<start_singlestriphit> & HStart2ClusterFinder::findcluster(std::vector<start_singlestriphit> & hitstoscan)
{
    std::vector<Bool_t> clustermask;
    clustermask.resize(hitstoscan.size(),0);
    std::vector<uint> clustersizes;

    for(UInt_t k = 0; k < hitstoscan.size(); k ++) {
	if(clustermask[k] == 1) { continue; }

	clustermask[k] = 1;
	std::vector<start_singlestriphit> timecluster;
	timecluster.resize(0);
	std::vector<start_singlestriphit> cluster4d;
	start_singlestriphit startstrip = hitstoscan[k];
	timecluster.push_back(startstrip);
	Int_t startstripnr = startstrip.getModID() * 100 + startstrip.getChannelID();
	std::vector<UInt_t> timeclusterindices;
	timeclusterindices.resize(0);
	Int_t currclusterwidth = 1;

	for(UInt_t i = k + 1; i < hitstoscan.size(); i ++) {
	    if(abs(hitstoscan[i].getTOA() - startstrip.getTOA()) < fToAthresh && clustermask[i] == 0) {
		timecluster       .push_back(hitstoscan[i]);
		timeclusterindices.push_back(i);
	    }
	}
	cluster4d.push_back(hitstoscan[k]);

	// form cluster in space
	assert(timeclusterindices.size() == timecluster.size() - 1);
	Bool_t consecutivecluster  = kTRUE;
	Bool_t abovemaxclustersize = kFALSE;
	std::vector<Int_t> clusterchannelnr;
	clusterchannelnr.push_back(startstripnr);
	for(UInt_t l = 0; l < timeclusterindices.size(); l ++) {
	    Int_t currstripnr = timecluster[l + 1].getModID() * 100 + timecluster[l + 1].getChannelID();
	    if(startstrip.getModID() == timecluster[l+1].getModID() && abs(currstripnr - startstripnr) < fStriprange) {
		currclusterwidth += 1;
		cluster4d.push_back(timecluster[l + 1]);
		clustermask[timeclusterindices[l]] = 1;
		clusterchannelnr.push_back(currstripnr);
	    }
	}
	std::sort(clusterchannelnr.begin(),clusterchannelnr.end());

	//discard clusters with non-consecutive strips
	for(UInt_t i = 1; i < clusterchannelnr.size(); i ++) {
	    if(clusterchannelnr[i] - clusterchannelnr[i - 1] != 1) { consecutivecluster = kFALSE; }
	}

	// apply clustercut
	if(cluster4d.size() > fMaxClustersize) { abovemaxclustersize = kTRUE; }

	if(consecutivecluster && !abovemaxclustersize) {

	    start_cluster currcluster;
	    currcluster.setclusterstrips(cluster4d);
	    fclustermaxstrips.push_back(currcluster.getmaxstripfromcluster());
	    fclustersizes    .push_back(currcluster.getClustersize());
	}
    }

    return fclustermaxstrips;
}
//###############################################################