//*-- Created : 03/08/2012 by S.Spataro
//_HADES_CLASS_DESCRIPTION
//////////////////////////////////////////////////////////////////////////////
//
//  HParticleStart2HitF
//
//  This class finds the correct start hit
//  Th correctiom method is store in HStart2Hit::gerCorrFlag() :
//      -1 : no corr, 0 : mod0, 1 : mod1. 2: both mods
//
///////////////////////////////////////////////////////////////////////////////

#include "hades.h"
#include "hcategory.h"
#include "hevent.h"
#include "hrecevent.h"
#include "hruntimedb.h"
#include "hspectrometer.h"
#include "htime.h"
#include "hbeamtime.h"
#include "hdatasource.h"
#include "hdetector.h"
#include "hitofdetector.h"

#include "hitofcal.h"
#include "hitofraw.h"
#include "hitofcalibrater.h"
#include "hitofcalpar.h"

#include "hparticlecand.h"
#include "hparticlecandsim.h"
#include "hstart2cal.h"
#include "hstart2hit.h"
#include "htofhit.h"
#include "htofcluster.h"
#include "hrpchit.h"
#include "hrpccluster.h"
#include "htofhitsim.h"
#include "htofclustersim.h"
#include "hrpchitsim.h"
#include "hrpcclustersim.h"
#include "hparticlestart2hitf.h"
#include "tofdef.h"
#include "rpcdef.h"
#include "hstartdef.h"
#include "hitofdef.h"
#include "hstart2clusterfinder.h"

#include "TMath.h"

#include <iomanip>
#include <iostream>

ClassImp(HParticleStart2HitF)

    HParticleStart2HitF::HParticleStart2HitF(void)
    : HReconstructor()
{
    // default constructor

    fitofRawCat  = NULL;
    fcandCat     = NULL;
    fCatTof      = NULL;
    fCatTofClu   = NULL;
    fCatRpc      = NULL;
    fCatRpcClu   = NULL;
    fCatStartCal = NULL;
    fCatStartHit = NULL;
    fpiTofCalPar = NULL;
    
    bDebug            = kFALSE;
    startTimeOriginal = 0.;
    startFlag         =-1;
    fVersion = 0;
    setTimeCut(15.);
}

HParticleStart2HitF::HParticleStart2HitF(const Text_t *name, const Text_t *title, Bool_t skip)
: HReconstructor(name, title)
{
    // constructor
    fitofRawCat  = NULL;
    fcandCat     = NULL;
    fpiTofCalPar = NULL;

    fCatTof      = NULL;
    fCatTofClu   = NULL;
    fCatRpc      = NULL;
    fCatRpcClu   = NULL;
    fCatStartCal = NULL;
    fCatStartHit = NULL;
    bDebug            = kFALSE;
    startTimeOriginal = 0.;
    startFlag         =-1;
    fVersion = 0;
    setTimeCut(5.,15.);
    memset(candCounter,0,sizeof(candCounter));
}

Bool_t HParticleStart2HitF::init(void)
{
    // init function

    if(!gHades) {
	Error("init()","No Hades object found!");
        return kFALSE;
    }

    //---------------------------------------------------------
    // switch the different beamtime versions

    if(!HBeamTime::checkBeamTimeID(gHades->getBeamTimeID(),"HParticleStart2HitF::init()",kTRUE)){ // HADES::kUnknownBeam rejected
	return kFALSE;
    }

    HBeamTime::checkBeamTimeSettingByRunId("HParticleStart2HitF::init()");


    Int_t year;
    HBeamTime::getYearFromRunID(year);

    if(year == 2022) {
	fVersion = 1;
        Info("init()", "Input file is detected from 2022. fVersion switch is set == 1 (LGAD cluster finder).");

    }

    if(year > 2022) { // needs future fixes! HBeamTime
	fVersion = 1;
        Info("init()", "Input file is detected after 2022. fVersion switch is set = 1 (LGAD cluster finder). If this is not desired, the code needs update");
    }
    //---------------------------------------------------------

    fCatTof = gHades->getCurrentEvent()->getCategory(catTofHit);
    if (NULL == fCatTof) {
	Warning("init", "HTofHit category not available!");
    }

    fCatTofClu = gHades->getCurrentEvent()->getCategory(catTofCluster);
    if (NULL == fCatTofClu) {
	Warning("init", "HTofCluster category not available!");
    }

    fCatRpc = gHades->getCurrentEvent()->getCategory(catRpcHit);
    if (NULL == fCatRpc) {
	Warning("init", "HRpcHit category not available!");
    }
    
    fitofRawCat  = gHades->getCurrentEvent()->getCategory(catiTofRaw);
    if (NULL == fitofRawCat) {
	Warning("init", "HiTofRawCat category not available!");
    }
    
    fcandCat = gHades->getCurrentEvent()->getCategory(catParticleCand);
    if (NULL == fcandCat) {
 	Warning("init", "HParticleCand category not available!");
    }
   
    fCatRpcClu = gHades->getCurrentEvent()->getCategory(catRpcCluster);
    if (NULL == fCatRpcClu) {
	Warning("init", "HRpcCluster category not available!");
    }

    fCatStartCal = gHades->getCurrentEvent()->getCategory(catStart2Cal);
    if (NULL == fCatStartCal) {
	Warning("init", "HStart2Cal category not available!");
    }

    fCatStartHit = gHades->getCurrentEvent()->getCategory(catStart2Hit);
    if (NULL == fCatStartHit) {
	Warning("init", "HStart2Hit category not available!");
    }

    fpiTofCalPar = (HiTofCalPar *)gHades->getRuntimeDb()->getContainer("iTofCalPar");

    if (bDebug)
    {
	out = new TFile("./start_output.root","RECREATE");
	nt1 = new TNtuple("nt1","nt1","evt:sys:tof:dist:tofn");
	nt2 = new TNtuple("nt2","nt2","evt:size0:size1:tof1:tof0:tof:n1:n0:n");
	nt3 = new TNtuple("nt3","nt3","evt:tof:start0:start1:dist0:dist1:startold:startnew:flag");
    }

    return kTRUE;
}

Int_t HParticleStart2HitF::execute(void)
{
    // makes the hits and fills the HStartHit category


    if(!fCatStartHit || !fCatStartCal) return 0;

    //--------------------------------------------------------------------
    // find the minimum tof from tof+rpc (already corrected for start hit)
    Float_t minTime = 0;
    if     (fVersion == 0) minTime = findMinimumTime(3);
    else if(fVersion == 1) minTime = findMinimumTime(1);
    else {
	Error("execute()","Unexpected fVersion = %i !",fVersion);
	exit(1);
    }

    if (minTime < 0) return 0; // if no good minimum tof value then skip the rest
    //--------------------------------------------------------------------

    //--------------------------------------------------------------------
    // find the start time from HStart2Hit
    startTimeOriginal = findOriginalStartTime();  // needs startHit
    if (startTimeOriginal == -1) return 0 ;
    //--------------------------------------------------------------------
    
    //--------------------------------------------------------------------
    // find new start time, "minTime + startTimeOriginal" is the time
    // before the start correction from start hit
    firstCluster.clear();
    secondCluster.clear();
    fSecondTime = -10000;

    Float_t startTime = 0;
    
    HStart2Hit* start2Hit = dynamic_cast<HStart2Hit*>(fCatStartHit->getObject(0));
    
    if (fVersion == 0) {
	startTime = findStartTime(minTime + startTimeOriginal);    // needs startCal
    } else if(fVersion == 1) {
	if (start2Hit->getIteration() == 1){
	    startTime = findStartTime_LGAD(minTime, startTimeOriginal);   // needs startCal
	} else if (start2Hit->getIteration() == 2 && fitofRawCat && fcandCat){
	    startFlag = start2Hit->getCorrFlag();
	    if(fitofRawCat && fcandCat){
		startTime = findStartTime_LGAD_itof_ParticleCand(start2Hit->getTime(),start2Hit->getCorrFlag());  // needs startCal,itofRaw,ParticleCand
	    } else  startTime  = startTimeOriginal;
	}
    } else {
	Error("execute()","Unexpected fVersion = %i !",fVersion);
	exit(1);
    }
    //--------------------------------------------------------------------


    //--------------------------------------------------------------------
    if (start2Hit->getIteration() == 2){
        correctParticleCand(startTime - startTimeOriginal);
    }
    
    correctTof(startTime - startTimeOriginal);
    correctRpc(startTime - startTimeOriginal);
    correctStart(startTime - startTimeOriginal);
    //--------------------------------------------------------------------

    return 0;
}

Float_t HParticleStart2HitF::findMinimumTime(Int_t numOfMetaHits)
{
    // Function which loops inside tof and rpc hits, find the 3 hits
    // with minimum time-of-flight, and returns the average value -> fastest signal

    HTofHit* pTofHit = NULL;
    HRpcHit* pRpcHit = NULL;

    vector<Float_t> tofList, rpcList, allList; // lists of time-of-flight

    //--------------------------------------------------------------------
    // loop in HTofHit category
    if(fCatTof) {
	for (Int_t tofEntry = 0; tofEntry < fCatTof->getEntries(); ++tofEntry)
	{
	    if (NULL == (pTofHit = static_cast<HTofHit*>(fCatTof->getObject(tofEntry)))) {
		Error("findMinimumTime", "Pointer to HTofHit object == NULL, returning");
		return -1;
	    }

	    //--------------------------------------------------------------------
            // skipp the embedded sim tracks
	    if(gHades->getEmbeddingMode()>0){
		HTofHitSim* pTofHitSim = 0;
		pTofHitSim = dynamic_cast<HTofHitSim*>(pTofHit);
		if(pTofHitSim){
                    if(pTofHitSim->getNTrack1() > 0 || pTofHitSim->getNTrack2() > 0) continue;
		}
	    }
	    //--------------------------------------------------------------------


	    Float_t tofTof = pTofHit->getTof();
	    if ( (tofTof < 0) || (tofTof > 300) ) continue; // cuts too low or high tof values

	    Float_t distTof = 0;
	    pTofHit->getDistance(distTof);
	    if (distTof > 3000) continue; // cuts not feasible distances

	    tofList.push_back(2100. * tofTof / distTof);  // fill with normalized tof to 2100 mm
	    allList.push_back(2100. * tofTof / distTof);  // fill with normalized tof to 2100 mm
	    if (bDebug) nt1->Fill(gHades->getEventCounter(), 1, tofTof, distTof, 2100. * tofTof / distTof);
	} // end of HTofHit loop
    }
    //--------------------------------------------------------------------

    //--------------------------------------------------------------------
    // loop in RpcHit category
    Int_t tracksL[10];
    Int_t tracksR[10];

    if(fCatRpc){
	for (Int_t rpcEntry = 0; rpcEntry < fCatRpc->getEntries(); ++rpcEntry)
	{
	    if (NULL == (pRpcHit = static_cast<HRpcHit*>(fCatRpc->getObject(rpcEntry)))) {
		Error("execute", "Pointer to HRpcHit object == NULL, returning");
		return -1;
	    }

	    //--------------------------------------------------------------------
            // skipp the embedded sim tracks
	    if(gHades->getEmbeddingMode()>0){
		HRpcHitSim* pRpcHitSim = 0;
		pRpcHitSim = dynamic_cast<HRpcHitSim*>(pRpcHit);
		if(pRpcHitSim){
		    pRpcHitSim->getTrackLArray(tracksL);
                    pRpcHitSim->getTrackRArray(tracksR);
                    Bool_t simTrack=kFALSE;
		    for(Int_t i = 0; i < 10; i++){
			if(tracksL[i] > 0 || tracksR[i] > 0) { simTrack = kTRUE; break;}
		    }
                    if(simTrack) continue;
		}
	    }
	    //--------------------------------------------------------------------

	    Float_t tofRpc = pRpcHit->getTof();
	    if ( (tofRpc < 0) || (tofRpc > 300) ) continue; // cuts too low or high tof values

	    Float_t xRpc = 0., yRpc = 0., zRpc = 0.;
	    pRpcHit->getXYZLab(xRpc, yRpc, zRpc);
	    Float_t distRpc = TMath::Sqrt(xRpc * xRpc + yRpc * yRpc + zRpc * zRpc);
	    if (distRpc > 3000) continue; // cuts not feasible distances

	    rpcList.push_back(2100. * tofRpc / distRpc);  // fill with normalized tof to 2100 mm;
	    allList.push_back(2100. * tofRpc / distRpc);  // fill with normalized tof to 2100 mm;
	    if (bDebug) nt1->Fill(gHades->getEventCounter(), 0, tofRpc, distRpc, 2100. * tofRpc / distRpc);
	} // end of HRpcHit loop
    }
    //--------------------------------------------------------------------

    if (allList.size() == 0) return -1; // if no good time, returns a negative value

    //--------------------------------------------------------------------
    // sort tof lists from minimum to maximum value
    std::sort(tofList.begin(), tofList.end());
    std::sort(rpcList.begin(), rpcList.end());
    std::sort(allList.begin(), allList.end());
    //--------------------------------------------------------------------


    //--------------------------------------------------------------------
    // select the three tof values with minimum time,
    // and calculate the average, for tof rpc and both
    Float_t tofSum   = 0, rpcSum   = 0, allSum   = 0;
    Float_t tofMean  = 0, rpcMean  = 0, allMean  = 0;
    Int_t   tofCount = 0, rpcCount = 0, allCount = 0;
    for (vector<Float_t>::iterator it = tofList.begin(); it != tofList.end(); ++it) {
	if (tofCount < numOfMetaHits) {
	    tofSum += *it;
	    tofCount++;
	}
    }
    if (tofCount > 0) tofMean = tofSum / tofCount;

    for (vector<Float_t>::iterator it = rpcList.begin(); it != rpcList.end(); ++it) {
	if (rpcCount < numOfMetaHits) {
	    rpcSum += *it;
	    rpcCount++;
	}
    }
    if (rpcCount > 0) rpcMean = rpcSum/rpcCount;

    for (vector<Float_t>::iterator it = allList.begin(); it != allList.end(); ++it) {
	if (allCount < numOfMetaHits) {
	    allSum += *it;
	    allCount++;
	}
    }
    if (allCount > 0) allMean = allSum / allCount;
    //--------------------------------------------------------------------

    if (bDebug) nt2->Fill(gHades->getEventCounter(), rpcList.size(), tofList.size(), tofMean, rpcMean, allMean, tofCount, rpcCount, allCount);

    return allMean; // returns the minimum time from tof+rpc
}

Float_t HParticleStart2HitF::findOriginalStartTime(void)
{
    // function to retrieve the start time from HStart2Hit category
    
    Float_t start_time = -1;
    //--------------------------------------------------------------------
    if(fCatStartHit)
    {
        HStart2Hit* start2Hit = dynamic_cast<HStart2Hit*>(fCatStartHit->getObject(0));
        if (!start2Hit) return -1; 
	    start_time = start2Hit->getTime();
    }
    //--------------------------------------------------------------------
    return start_time;
}

Float_t HParticleStart2HitF::findStartTime_LGAD(Float_t minTof, Float_t oldStartHitTime)
{
    // Find the correct start time, where time-start closer to 7ns
    
    Float_t startTimeDiff[2] = { 300., 300.};
    Float_t startTime    [2] = { 300., 300. };
    startFlag = -1;

    std::vector<start_singlestriphit> allstriphitsperentry;
    allstriphitsperentry.clear();
    if(fCatStartCal){
	HStart2Cal* pStartCal = NULL;
	for (Int_t startEntry = 0; startEntry < fCatStartCal->getEntries(); ++startEntry){
	    // loop on HStartCal
	    if (NULL == (pStartCal = static_cast<HStart2Cal*>(fCatStartCal->getObject(startEntry)))) {
		Error("findStartTime", "Pointer to HStartCal object == NULL, returning");
		return 0;
	    }
	    Int_t mod       = pStartCal->getModule();
	    Int_t currStrip = pStartCal->getStrip();
	    if (mod > 1) continue; // skip VETO
	    for (Int_t nStart = 0; nStart < pStartCal->getMultiplicity(); nStart++){ // loop on the start time
		Float_t time  = pStartCal->getTime (nStart + 1);
		Float_t width = pStartCal->getWidth(nStart + 1);
		
		start_singlestriphit currsinglestriphit(mod,currStrip,nStart, time ,width);
		allstriphitsperentry.push_back(currsinglestriphit);
	    }
	    
	}
    }
    // cluster finding
    HStart2ClusterFinder clfinder;
    clfinder.setclusterthresholds(2,5, 3); //ToArange,striprange,maxclustersizeinstrips
    std::vector<start_singlestriphit> localMaxStrips = clfinder.findcluster(allstriphitsperentry);
    
    
    
    //--------------------------------------------------------------------
    
    if(localMaxStrips.size() != 0){
	for (UInt_t i = 0; i < localMaxStrips.size(); i ++){
	    Int_t mod = localMaxStrips[i].getModID();
	    Float_t time = localMaxStrips[i].getTOA();
	    if (( time < -20) || ( time > 50) ) continue;
	    
	    Float_t distanceTarget_Meta = 2100.;
	    Float_t distanceTarget_iTOF = 650.;// TODO estimate real distance between target and iTOF
	    
	    Float_t metaHitTime            = minTof + oldStartHitTime ;
	    Float_t expectedT0_numerator   = - ((metaHitTime - 7.) * distanceTarget_iTOF - oldStartHitTime * distanceTarget_Meta);
	    Float_t expectedT0_denominator = (distanceTarget_Meta - distanceTarget_iTOF);
	    
	    Float_t expectedT0 = expectedT0_numerator / expectedT0_denominator;
	    
	    
	    if (TMath::Abs(expectedT0 - time ) < TMath::Abs(startTimeDiff[mod]) && TMath::Abs(expectedT0 - time ) < 3 ){
		startTimeDiff[mod] = expectedT0 - time;
		startTime    [mod] = time;
	    }
	}
    }
    

    //--------------------------------------------------------------------
    Float_t newStart = oldStartHitTime;
    if ( (startTime[0] != 300) && (startTime[1] != 300) ) // if there is good start time for both modules
    {
	// if they are within 0.5ns then calculates mean value
	if (TMath::Abs(startTime[0] - startTime[1]) < 0.5)
	{
	    newStart = 0.5 * (startTime[0] + startTime[1]);
	    startFlag = 2;
	} else {
	    // if the 2 start times are not correlated, take the closer to 7 ns
	    if (TMath::Abs(startTimeDiff[0]) < TMath::Abs(startTimeDiff[1]))
	    {
		newStart = startTime[0];
		startFlag = 0;
	    }
	    else
	    {
		newStart = startTime[1];
		startFlag = 1;
	    }
	}
    } // end of both modules
    else
    {
	// if only one start time is present, select such module for time
	if (startTime[0] != 300)
	{
	    newStart = startTime[0];
	    startFlag = 0;
	}
	if (startTime[1] != 300)
	{
	    newStart = startTime[1];
	    startFlag = 1;
	}
    }
    //--------------------------------------------------------------------


    // if no good time is present in no start modules, the start time remains "0"

    if (bDebug) nt3->Fill(gHades->getEventCounter(), minTof, startTime[0], startTime[1], startTimeDiff[0], startTimeDiff[1], startTimeOriginal, newStart, startFlag);

    return newStart;
}

Int_t HParticleStart2HitF::findiTOFIndex_Raw(Int_t sector, Double_t phi, Float_t oldStartHitTime, Float_t originalToF,Float_t theta )
{
    Int_t iTOFIndex = -1000;
    Int_t numberOfEntriesInSamePhi = 0;
    Float_t iTOF_phiDistribution[18] = {75,90,105,135,150,165,195,210,225,255,270,285,315,330,345,15,30,45};
    for (Int_t i = 0; i < fitofRawCat->getEntries();i++){
	HiTofRaw*  itofRaw = (HiTofRaw *)fitofRawCat ->getObject(i);
	if (!itofRaw){ continue; }
	Int_t itofSector = itofRaw->getSector();
	Int_t itofPad    = itofRaw->getCell();

	Float_t phiDifference = phi - iTOF_phiDistribution[itofSector*3+itofPad];
	if (phiDifference >= 180.) phiDifference -=360.;
	if (phiDifference <=-180.) phiDifference +=360.;
	if (fabs(phiDifference) < 8 && itofRaw->getSector() == sector){
	    if (numberOfEntriesInSamePhi == 0){
		numberOfEntriesInSamePhi++;
		iTOFIndex = i;
	    }
	    else if (numberOfEntriesInSamePhi == 1) iTOFIndex = -1000;
	}
    }
    return iTOFIndex;
}
    
    
    
void HParticleStart2HitF::correctParticleCand(Float_t t0_corretion)
{

    // is a two step procedure; first recalculate all betas+mass,
    // then use them to modify the category values. Othwerwise,
    // the modified beta value would modify the output of the
    // loop if done in just one single step!
    if(!fcandCat) return;


    HParticleCand* pCand=0;
    if(fcandCat)
    {
	for(Int_t j = 0; j < fcandCat->getEntries(); j++){
	    pCand = dynamic_cast<HParticleCand*>(fcandCat->getObject(j));
	    if(!pCand) continue;
	    Float_t c = 299.792458;
	    Int_t systemUsed =-1;
	    systemUsed=pCand->getSystemUsed();
	    Int_t rpcClusterIndex =-1;
	    Int_t tofClusterIndex = -1;
	    Float_t oldTofMeta = -10.;

	    HTofCluster* pTofCluster = NULL;
	    HRpcCluster* pRpcCluster = NULL;
	    if (systemUsed == 0 ){
		rpcClusterIndex=pCand->getRpcInd();
		if (rpcClusterIndex != -1){
		    pRpcCluster = dynamic_cast<HRpcCluster*>(fCatRpcClu->getObject(rpcClusterIndex));
		    oldTofMeta  = pRpcCluster->getTof();
		}
	    }
	    else if (systemUsed == 1){
		tofClusterIndex = pCand->getTofClstInd();
		if (tofClusterIndex != -1){
		    pTofCluster = dynamic_cast<HTofCluster*>(fCatTofClu->getObject(tofClusterIndex));
		    oldTofMeta  = pTofCluster->getTof();
		}
	    }
	    else if ( tofClusterIndex == -1 && rpcClusterIndex && oldTofMeta == -10) continue;

	    // TODO check edge conditions of getTof / getBeta
	    Float_t newBeta = pCand->getDistanceToMetaHit()/(oldTofMeta - t0_corretion)/c;
	    if ((oldTofMeta - t0_corretion) > 0) pCand->setBeta(newBeta);
	    else  pCand->setBeta(-1);
	    Float_t mom = pCand->getMomentum();

	    Float_t mass2 = mom*mom/newBeta/newBeta*(1.- (newBeta*newBeta)  );
	    if (mass2>0) pCand->setMass2(mass2);
	    else mass2 = -10000;

	}
    }
}

    
void HParticleStart2HitF::countParticlesInItofPaddles()
{

    if (NULL == fcandCat) {
	return;
    }

    Float_t iTOF_phiDistribution[18]= {75,90,105,135,150,165,195,210,225,255,270,285,315,330,345,15,30,45};
    for (Int_t i = 0; i < fcandCat->getEntries(); i++){
	HParticleCand* cand = dynamic_cast <HParticleCand*> (fcandCat->getObject(i));
	if (!cand) continue;
	for (Int_t j = 0;j < 18; j++){ // select paddle index
	    Float_t phiDifference = cand->getPhi() - iTOF_phiDistribution[j];
	    if (phiDifference >= 180.)phiDifference -= 360.;
	    if (phiDifference <=-180.)phiDifference += 360.;
	    if (fabs(phiDifference) <8 ){ candCounter[j]++; }
	}
    }
}

Bool_t customComparerForSorting(const std::pair<Float_t, Float_t> &p1, const std::pair<Float_t, Float_t> &p2)
{
  return p1.first < p2.first ;
  
}    
    
Float_t HParticleStart2HitF::findStartTime_LGAD_itof_ParticleCand(Float_t oldStartHitTime, Int_t oldStartFlag)
{
    // Find the correct start time, where time-start closer to 7ns

    if (NULL == fcandCat) return oldStartHitTime;

    Float_t startTimeDiff[2] = { 300., 300.};
    Float_t startTime    [2] = { 300., 300. };

    Float_t itofDistanceParabola[6][3]= {   {772.156, -8.98532, 0.100889},
    {773.271, -8.99527, 0.101},
    {776.583, -9.00908, 0.101494},
    {780.18, -9.00421,  0.102202},
    {779.355, -9.00543, 0.102039},
    {778.824, -9.02441, 0.101772} };

    std::vector<start_singlestriphit> allstriphitsperentry;
    allstriphitsperentry.clear();
    if(fCatStartCal){
	HStart2Cal* pStartCal = NULL;

	for (Int_t startEntry = 0; startEntry < fCatStartCal->getEntries(); ++startEntry){
	    // loop on HStartCal
	    if (NULL == (pStartCal = static_cast<HStart2Cal*>(fCatStartCal->getObject(startEntry)))) {
		Error("findStartTime_LGAD_itof_ParticleCand", "Pointer to HStartCal object == NULL, returning");
		return 0;
	    }
	    Int_t mod       = pStartCal->getModule();
	    Int_t currStrip = pStartCal->getStrip();
	    if (mod > 1) continue; // skip VETO
	    for (Int_t nStart = 0; nStart < pStartCal->getMultiplicity(); nStart++){ // loop on the start time
		Float_t time  = pStartCal->getTime (nStart + 1);
		Float_t width = pStartCal->getWidth(nStart + 1);

		start_singlestriphit currsinglestriphit(mod,currStrip,nStart, time ,width);
		allstriphitsperentry.push_back(currsinglestriphit);
	    }

	}
	if (allstriphitsperentry.size()==0){
	}
    }

    // cluster finding
    HStart2ClusterFinder clfinder;
    clfinder.setclusterthresholds(2,5, 3); //ToArange,striprange,maxclustersizeinstrips
    std::vector<start_singlestriphit> localMaxStrips = clfinder.findcluster(allstriphitsperentry);

    //--------------------------------------------------------------------

    Float_t sumExpectedT0 = -1000;
    Int_t counterExpectedT0 = 0;
    HParticleCand *cand = NULL;
    Int_t sizeCandCat = fcandCat->getEntries();
    Float_t theta = -1000;
    Float_t fastestTof = 10000;
    Float_t bestchi2 = 1000000.;
    Float_t expectedT0_bestchi2 = -1000;
    Float_t earliest_expectedT0 = 10000;

    memset(candCounter,0,sizeof(candCounter));

    countParticlesInItofPaddles();
    Int_t counterCandSurvivingSelection = 0;

    std::vector<std::pair<Float_t, Float_t> > vector_chi2_and_expectedT0;

    for (Int_t i = 0; i < sizeCandCat; i++){
	cand = dynamic_cast <HParticleCand*> (fcandCat->getObject(i));
	if (!cand) continue;
	if (  ! ((cand->getChi2() > 0) && (cand->getChi2() < 1000))  )continue;

	if (cand->getDistanceToMetaHitOrg() == 0) continue;
	if (cand->getTof() == -1)                 continue;
	counterCandSurvivingSelection++;
	if (cand->getTheta() == -1000)            continue;
	theta = cand->getTheta();

	Int_t iTOFIndex = findiTOFIndex_Raw (cand->getSector(),cand->getPhi(), oldStartHitTime , cand->getTof(),theta) ;

	if (iTOFIndex == -1000) continue;
	HiTofRaw* itofRaw = (HiTofRaw *)fitofRawCat ->getObject(iTOFIndex);
	if (!itofRaw)           continue;

	Float_t fast_time = 10000.;
	Float_t sumTime = 0;
	Int_t timeCounter = 0;

	for (Int_t j = 0; j < 12; j++){// loop over PMTs
	    if (itofRaw->getTime(j) <= -100000) continue;
	    Float_t calRawTime    = fpiTofCalPar->calibrateTime_walkcorr           (itofRaw->getSector(), itofRaw->getModule(), itofRaw->getCell(), j, itofRaw->getTime(j),itofRaw->getWidth(j));
	    Float_t newCalRawTime = fpiTofCalPar->correctSignalPropagation_walkcorr(itofRaw->getSector(), itofRaw->getModule(), itofRaw->getCell(), j, calRawTime,theta);
	    if (newCalRawTime<fast_time && fast_time > -20 )  fast_time =   newCalRawTime;
	    sumTime += newCalRawTime;
	    timeCounter++;
	}

	Float_t metaHitTime  = cand->getTof() + oldStartHitTime ;
	Float_t itof_hitTime = fast_time ;
	Float_t distanceTarget_iTOF  =  itofDistanceParabola[itofRaw->getSector()][0] + (itofDistanceParabola[itofRaw->getSector()][1]*cand->getTheta())
	    + (itofDistanceParabola[itofRaw->getSector()][2]*cand->getTheta()*cand->getTheta() );

	Float_t expectedT0_numerator   = -((metaHitTime) * (distanceTarget_iTOF) - itof_hitTime * cand->getDistanceToMetaHit());
	Float_t expectedT0_denominator = ((cand->getDistanceToMetaHitOrg()) - (distanceTarget_iTOF));
	Float_t expectedT0 = expectedT0_numerator / expectedT0_denominator;

	std::pair<Float_t,Float_t> chi2_and_expectedT0;

	chi2_and_expectedT0.first  = cand->getChi2();
	chi2_and_expectedT0.second = expectedT0;
	vector_chi2_and_expectedT0.push_back(chi2_and_expectedT0);

	if (cand->getTof() != -1){
	    if (cand->getTof() < fastestTof){
		fastestTof = cand->getTof();
	    }
	    if (cand->getChi2() < bestchi2){
		bestchi2 = cand->getChi2();
		expectedT0_bestchi2 = expectedT0;
	    }

	    if (expectedT0 < earliest_expectedT0) earliest_expectedT0 = expectedT0;
	}

	if (sumExpectedT0 == -1000) sumExpectedT0 = expectedT0;
	else sumExpectedT0 += expectedT0;
	counterExpectedT0++;
    }

    std::stable_sort(vector_chi2_and_expectedT0.begin(),vector_chi2_and_expectedT0.end(), customComparerForSorting) ;

    Float_t expectedT0 = -1000;
    Float_t timeCut = 3.;
    Bool_t foundT0 = kFALSE;
    for (UInt_t j = 0; j < vector_chi2_and_expectedT0.size(); j++){
	if (foundT0)continue;
	expectedT0 = vector_chi2_and_expectedT0[j].second;
	if(localMaxStrips.size() != 0 && expectedT0 != -1000){
	    for (UInt_t i = 0; i < localMaxStrips.size(); i ++){
		Int_t mod = localMaxStrips[i].getModID();
		Float_t time = localMaxStrips[i].getTOA();
		if (( time < -20) || ( time > 50) ) continue;

		if (TMath::Abs(expectedT0 - time ) < TMath::Abs(startTimeDiff[mod]) && TMath::Abs(expectedT0 - time ) < (timeCut +20)  ){
		    startTimeDiff[mod] = expectedT0 - time;
		    startTime    [mod] = time;
		    if (fabs(startTimeDiff[mod])<timeCut){
			foundT0=kTRUE;
		    }
		}
	    }
	}
    }

    //--------------------------------------------------------------------
    Float_t newStart = oldStartHitTime;
    if ( (startTime[0] < 300) && (startTime[1] != 300) ) // if there is good start time for both modules
    {
	// if they are within 0.5ns then calculates mean value
	if (TMath::Abs(startTime[0] - startTime[1]) < 0.5  && ( TMath::Abs(expectedT0-startTime[0]) < timeCut || TMath::Abs(expectedT0-startTime[1]) < timeCut      ) )
	{
	    newStart = 0.5 * (startTime[0] + startTime[1]);
	    if (newStart == oldStartHitTime) startFlag = 35;
	    else startFlag = 30+oldStartFlag+1;
	} else {
	    // if the 2 start times are not correlated, take the closer to 7 ns
	    if (TMath::Abs(startTimeDiff[0]) < TMath::Abs(startTimeDiff[1])  && ( TMath::Abs(startTime[0] - expectedT0) < timeCut  ) )
	    {
		newStart = startTime[0];
		if (newStart == oldStartHitTime) startFlag = 15;
		else startFlag = 10+oldStartFlag+1;
	    }
	    else if (TMath::Abs(startTime[1] - expectedT0) < timeCut)
	    {
		newStart = startTime[1];
		if (newStart == oldStartHitTime) startFlag = 25;
		else startFlag = 20+oldStartFlag+1;
	    }
	}
    } // end of both modules
    else
    {
	// if only one start time is present, select such module for time
	if (startTime[0] != 300 &&  TMath::Abs(expectedT0-startTime[0]) < timeCut )
	{
	    newStart = startTime[0];
	    if (newStart == oldStartHitTime) startFlag = 15;
	    else startFlag = 10+oldStartFlag+1;
	}
	if (startTime[1] != 300 &&  TMath::Abs(expectedT0-startTime[1]) < timeCut )
	{
	    newStart = startTime[1];
	    if (newStart == oldStartHitTime) startFlag = 25;
	    else startFlag = 20+oldStartFlag+1;
	}
    }
    //--------------------------------------------------------------------

    if (startFlag == -1) newStart = expectedT0_bestchi2;

    return newStart;
}

Float_t HParticleStart2HitF::findStartTime(Float_t minTof)
{
    // Find the correct start time, where time-start closer to 7ns

    
    Float_t startTimeDiff[2] = { 300., 300.};
    Float_t startTime    [2] = { 300., 300. };

    Int_t strip[2]={-1,-1};
    //--------------------------------------------------------------------
    if(fCatStartCal)
    {
	HStart2Cal* pStartCal = NULL;
	for (Int_t startEntry = 0; startEntry < fCatStartCal->getEntries(); ++startEntry)
	{
	    // loop on HStartCal
	    if (NULL == (pStartCal = static_cast<HStart2Cal*>(fCatStartCal->getObject(startEntry)))) {
		Error("findStartTime", "Pointer to HStartCal object == NULL, returning");
		return -1;
	    }

	    Int_t mod = pStartCal->getModule();
	    if (mod > 1) continue; // skip VETO

	    for (Int_t nStart = 0; nStart < pStartCal->getMaxMultiplicity(); nStart++) // loop on the start time
	    {
		Float_t time = pStartCal->getTime(nStart + 1);
		//if ((time<-10) || (time>230) ) continue;  // use only [-10,230] ns as confidential range
		if ( ((minTof-time) < timeCutL) || ((minTof-time) > timeCutR) ) continue; // consider start hits where tof is [5,15] ns

		// select start time (modulewise) where tof-start is closer to 7 ns.
		if (TMath::Abs(minTof - time - 7.) < TMath::Abs(startTimeDiff[mod]))
		{
		    startTimeDiff[mod] = minTof - time - 7.;
		    startTime    [mod] = time;
                    strip        [mod] = startEntry*1000+nStart+1;
		}
	    } // end of loop inside 4 start time

	} // end of start hit loop
    }
    //--------------------------------------------------------------------

    //--------------------------------------------------------------------
    Float_t newStart = 0;
    startFlag = -1;
    if ( (startTime[0] != 300) && (startTime[1] != 300) ) // if there is good start time for both modules
    {
	// if they are within 0.5ns then calculates mean value
	if (TMath::Abs(startTime[0] - startTime[1]) < 0.5)
	{
	    newStart = 0.5 * (startTime[0] + startTime[1]);
	    startFlag = 2;
	} else {
	    // if the 2 start times are not correlated, take the closer to 7 ns
	    if (TMath::Abs(startTimeDiff[0]) < TMath::Abs(startTimeDiff[1]))
	    {
		newStart = startTime[0];
		startFlag = 0;
	    }
	    else
	    {
		newStart = startTime[1];
		startFlag = 1;
	    }
	}
    } // end of both modules
    else
    {
	// if only one start time is present, select such module for time
	if (startTime[0] != 300)
	{
	    newStart = startTime[0];
	    startFlag = 0;
	}
	if (startTime[1] != 300)
	{
	    newStart = startTime[1];
	    startFlag = 1;
	}
    }
    //--------------------------------------------------------------------

    //--------------------------------------------------------------------
    // collect cluster infos
    if(startFlag>-1){ // any good found

        Float_t minDiffSecond = 0.5;

        // put the closest strip allways first in vector
        if(startFlag == 0 || startFlag == 2) firstCluster.push_back(strip[0]);
        if(startFlag == 1 || startFlag == 2) firstCluster.push_back(strip[1]);

        vector<Int_t> vtemp;   // temp vector for strips will be refilled later
	Int_t    nextHit =-1;
        Float_t timediff = 100000;   // local var
        Float_t newStart2= -100000;  // clostest scond time
	HStart2Cal* pStartCal = NULL;
	for (Int_t startEntry = 0; startEntry < fCatStartCal->getEntries(); ++startEntry)
	{
	    // loop on HStartCal
	    if (NULL == (pStartCal = static_cast<HStart2Cal*>(fCatStartCal->getObject(startEntry)))) {
		Error("findStartTime", "Pointer to HStartCal object == NULL, returning");
		return -1;
	    }

	    Int_t mod = pStartCal->getModule();
	    if (mod > 1) continue; // skip VETO

	    for (Int_t nStart = 0; nStart < pStartCal->getMaxMultiplicity(); nStart++) // loop on the 4 start time
	    {
		Float_t time = pStartCal->getTime(nStart + 1);
		//if ((time<-10) || (time>230) ) continue;  // use only [-10,230] ns as confidential range
		if ( ((minTof-time) < timeCutL) || ((minTof-time) > timeCutR) ) continue; // consider start hits where tof is [5,15] ns
                Int_t newstrip = startEntry*1000+nStart+1 ;
		// fill all strips of first cluster

                Float_t diff =  TMath::Abs( newStart - time );
		if (diff < 0.5)
		{
		    if(find(firstCluster.begin(),firstCluster.end(),newstrip)==firstCluster.end()){
			firstCluster.push_back(newstrip);
		    }

		} else {

		    if (diff > minDiffSecond)
		    {
			if(diff<timediff) {
			    timediff  = diff;
			    newStart2 = time;
			    nextHit   = newstrip;
			}
			vtemp.push_back(newstrip);
		    }
		}
	    } // end of loop inside 4 start time

	} // end of start hit loop

	if(vtemp.size()>0){  // there was something else
	    fSecondTime = newStart2;
	    secondCluster.push_back(nextHit);
	    for(UInt_t i=0;i<vtemp.size();i++){
		if(find(secondCluster.begin(),secondCluster.end(),vtemp[i])==secondCluster.end()){
		    secondCluster.push_back(vtemp[i]);
		}
	    }
	}

    }

    //--------------------------------------------------------------------

    // if no good time is present in no start modules, the start time remains "0"

    if (bDebug) nt3->Fill(gHades->getEventCounter(), minTof, startTime[0], startTime[1], startTimeDiff[0], startTimeDiff[1], startTimeOriginal, newStart, startFlag);

    return newStart;
}

void HParticleStart2HitF::correctTof(Float_t corrTime)
{
    // function which loops into HTofHit/HTofCluster categories and correct the tof
    // newtof = oldtof - corrTime

    //--------------------------------------------------------------------
    if(fCatTof){
	HTofHit* pTofHit = NULL;
	for (Int_t tofEntry = 0; tofEntry < fCatTof->getEntries(); ++tofEntry)
	{
	    // loop in HTofHit category
	    if (NULL == (pTofHit = static_cast<HTofHit*>(fCatTof->getObject(tofEntry)))) {
		Error("correctTof", "Pointer to HTofHit object == NULL, returning");
		return;
	    }

	    //--------------------------------------------------------------------
            // skipp the embedded sim tracks
	    if(gHades->getEmbeddingMode()>0){
		HTofHitSim* pTofHitSim = 0;
		pTofHitSim = dynamic_cast<HTofHitSim*>(pTofHit);
		if(pTofHitSim){
                    if(pTofHitSim->getNTrack1() > 0 || pTofHitSim->getNTrack2() > 0) continue;
		}
	    }
	    //--------------------------------------------------------------------

	    Float_t tofTof = pTofHit->getTof();
	    pTofHit->setTof(tofTof - corrTime);
	} // end of HTofHit loop
    }
    //--------------------------------------------------------------------

    //--------------------------------------------------------------------
    if(fCatTofClu){
	HTofCluster* pTofCluster = NULL;
	for (Int_t tofEntry = 0; tofEntry < fCatTofClu->getEntries(); ++tofEntry)
	{
	    // loop in HTofCluster category
	    if (NULL == (pTofCluster = static_cast<HTofCluster*>(fCatTofClu->getObject(tofEntry)))) {
		Error("findMinimumTime", "Pointer to HTofCluster object == NULL, returning");
		return;
	    }

	    //--------------------------------------------------------------------
            // skipp the embedded sim tracks
	    if(gHades->getEmbeddingMode()>0){
		HTofClusterSim* pTofClusterSim = 0;
		pTofClusterSim = dynamic_cast<HTofClusterSim*>(pTofCluster);
		if(pTofClusterSim){
                    if(pTofClusterSim->getNTrack1() > 0 || pTofClusterSim->getNTrack2() > 0) continue;
		}
	    }
	    //--------------------------------------------------------------------

	    Float_t tofTof = pTofCluster->getTof();
	    pTofCluster->setTof(tofTof - corrTime);
	} // end of HTofCluster loop
    }
    //--------------------------------------------------------------------

    return;
}

void HParticleStart2HitF::correctRpc(Float_t corrTime)
{
    // function which loops into HRpcHit/HRpcCluster categories and correct the tof
    // newtof = oldtof - corrTime

    //--------------------------------------------------------------------

    if(fCatRpc){
        Int_t tracksL[10];
        Int_t tracksR[10];

	HRpcHit* pRpcHit = NULL;
	for (Int_t tofEntry = 0; tofEntry < fCatRpc->getEntries(); ++tofEntry)
	{
	    // loop in HRpcHit category
	    if (NULL == (pRpcHit = static_cast<HRpcHit*>(fCatRpc->getObject(tofEntry)))) {
		Error("correctRpc", "Pointer to HRpcHit object == NULL, returning");
		return;
	    }
	    //--------------------------------------------------------------------
            // skipp the embedded sim tracks
	    if(gHades->getEmbeddingMode()>0){
		HRpcHitSim* pRpcHitSim = 0;
		pRpcHitSim = dynamic_cast<HRpcHitSim*>(pRpcHit);
		if(pRpcHitSim){
		    pRpcHitSim->getTrackLArray(tracksL);
                    pRpcHitSim->getTrackRArray(tracksR);
                    Bool_t simTrack=kFALSE;
		    for(Int_t i = 0; i < 10; i++){
			if(tracksL[i] > 0 || tracksR[i] > 0) { simTrack = kTRUE; break;}
		    }
                    if(simTrack) continue;
		}
	    }
	    //--------------------------------------------------------------------

	    Float_t tofRpc  = pRpcHit->getTof();
	    Float_t newTof  = tofRpc - corrTime;
	    Float_t acharge = pRpcHit->getCharge();
	    Float_t axmod   = pRpcHit->getXMod();
	    Float_t aymod   = pRpcHit->getYMod();
	    Float_t azmod   = pRpcHit->getZMod();

	    pRpcHit->setHit(newTof, acharge, axmod, aymod, azmod);
	} // end of HRpcHit loop
    }
    //--------------------------------------------------------------------

    //--------------------------------------------------------------------
    if(fCatRpcClu){
        Int_t tracks[4];
	HRpcCluster* pRpcCluster = NULL;
	for (Int_t tofEntry = 0; tofEntry < fCatRpcClu->getEntries(); ++tofEntry)
	{
	    // loop in HRpcCluster category
	    if (NULL == (pRpcCluster = static_cast<HRpcCluster*>(fCatRpcClu->getObject(tofEntry)))) {
		Error("findMinimumTime", "Pointer to HRpcCluster object == NULL, returning");
		return;
	    }

	    //--------------------------------------------------------------------
            // skipp the embedded sim tracks
	    if(gHades->getEmbeddingMode()>0){
		HRpcClusterSim* pRpcClusterSim = 0;
		pRpcClusterSim = dynamic_cast<HRpcClusterSim*>(pRpcCluster);
		if(pRpcClusterSim){
		    pRpcClusterSim->getTrackList(tracks);
                    Bool_t simTrack=kFALSE;
		    for(Int_t i = 0; i < 4; i++){
			if(tracks[i] > 0) { simTrack = kTRUE; break;}
		    }
                    if(simTrack) continue;
		}
	    }
	    //--------------------------------------------------------------------

	    Float_t tofRpc  = pRpcCluster->getTof();
	    Float_t newTof  = tofRpc - corrTime;
	    Float_t acharge = pRpcCluster->getCharge();
	    Float_t axmod   = pRpcCluster->getXMod();
	    Float_t aymod   = pRpcCluster->getYMod();
	    Float_t azmod   = pRpcCluster->getZMod();

	    pRpcCluster->setCluster(newTof, acharge, axmod, aymod, azmod);
	} // end of HRpcCluster loop
    }
    //--------------------------------------------------------------------

    return;
} // end of correctRpc function

void HParticleStart2HitF::correctStart(Float_t corrTime)
{
    // function which loops into HStart2Hit category and correct the tof
    // newtof = oldtof - corrTime

    //--------------------------------------------------------------------
    if(fCatStartHit){
	HStart2Hit* pStartHit = NULL;
	for (Int_t startEntry = 0; startEntry < fCatStartHit->getEntries(); ++startEntry)
	{
	    // loop in HStart2Hit category
	    if (NULL == (pStartHit = static_cast<HStart2Hit*>(fCatStartHit->getObject(startEntry)))) {
		Error("correctStart", "Pointer to HStart2Hit object == NULL, returning");
		return;
	    }

	    Float_t timeStart = 0, widthStart = 0;
	    pStartHit->getTimeAndWidth(timeStart, widthStart);
	    pStartHit->setTimeAndWidth(timeStart + corrTime, widthStart);
        
        if (startFlag < 5){
            pStartHit->setCorrFlag(startFlag);
            pStartHit->setIteration(2);
        }
        else if (startFlag >5 ){
            pStartHit->setIteration(3);
            pStartHit->setCorrFlag(startFlag/10 -1);
        }

	    //--------------------------------------------------------------------
	    // fill cluster infos
	    if(startFlag!=-1){
		pStartHit->resetClusterStrip(2);
		for(UInt_t i=0;i<firstCluster.size();i++){
		    pStartHit->setClusterStrip (0, firstCluster[i]);
		}
		for(UInt_t i=0;i<secondCluster.size();i++){
		    pStartHit->setClusterStrip (1, secondCluster[i]);
		}
		pStartHit->setSecondTime(fSecondTime);

		if(secondCluster.size()>0) pStartHit->setMultiplicity(2);
                else                       pStartHit->setMultiplicity(1);
	    }
	    //--------------------------------------------------------------------

	} // end of HStart2Hit loop
    }
    //--------------------------------------------------------------------

    return;
}

Bool_t HParticleStart2HitF::finalize(void)
{
    if (bDebug){
	out->cd();
	nt1->Write();
	nt2->Write();
	nt3->Write();
	out->Save();

	nt1->Delete();
	nt2->Delete();
	nt3->Delete();
	nt1 = 0;
	nt2 = 0;
	nt3 = 0;

	out->Close();
    }

    return kTRUE;
}