#include <iostream>
#include <cstdlib>
#include <cmath>
#include <bitset> 

#include <TGraph.h>
#include <TFile.h>

#include "TVector3.h"

#ifndef NOSKLIBRARIES
#include "skheadC.h"
#else
#define COUNT_PER_NSEC (1.92000)
#endif

#include "fQEventManager.h"

#include "spliTChan.h"
#include "spliTChanOutC.h"
#include "TRuntimeParameters.hxx"

extern"C" {
  
  // ${CERN}/2005/src/packlib/kernlib/kerngen/tcgen/sortzv.F
  // http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/m101/top.html
  void sortzv_(float *, int *, int *, int *, int *, int *);
  
}

int spliTChan::fMAXPMTZ = 334380;
int spliTChan::fNPMT = 11146; 
int spliTChan::fNMULT = 30;

using namespace std;

spliTChan::spliTChan(int aWriteHistos, int anPMT, int anMult ){
  
  writeHistos = aWriteHistos;

  fMAXPMTZ = anPMT*anMult;
  fNPMT = anPMT;
  fNMULT = anMult;

  // "Constant" parameters
  q_thresh = 0.3; 
  nwtr = 1.38;
  cwt = 30/nwtr;

  doTOFcorr = 1;

  //t_ind.reserve( fMAXPMTZ );
  //t_ind.resize( fMAXPMTZ );

  t_ind = new int [ fMAXPMTZ ];
  chrgz = new float [ fMAXPMTZ ];
  tHitz = new float [ fMAXPMTZ ];

  t = new float [ fMAXPMTZ ];
  q = new float [ fMAXPMTZ ];

  //amis
  ttc[1] = new float [ fMAXPMTZ ];
  ttc[0] = new float [ fMAXPMTZ ];

  for (int ic=0; ic<MAXCAND; ic++) chrg[ic] = new float [fMAXPMTZ];
  for (int ic=0; ic<MAXCAND; ic++) tHit[0][ic] = new float[fMAXPMTZ];
  for (int ic=0; ic<MAXCAND; ic++) tHit[1][ic] = new float[fMAXPMTZ];
  
  h_tiskz = 0;
  for (int icand=0; icand<MAXCAND; icand++) {
    h_tisk_corr0[icand] = 0;
    h_tisk_corr1[icand] = 0;
  }
  
  Cluster_tEnd = new float[MAXCAND];
  Cluster_tStart = new float[MAXCAND];
  //Cluster_nPeaks = new int *[MAXCAND];
  //for (int icluster=0; icluster<MAXCAND; icluster++)
  //Cluster_nPeaks[icluster] = new int[NPEAKFINDS];
  Cluster_Init();

  // Cluster search parameters //
  Cluster_nBinsScanWin = 7; // *10 ns
  Cluster_nHitThr = 9;
  Cluster_nHitThrFall = 9;
  Cluster_tBinSepThresh = 3;
  //Cluster_tExtendMin = 100;
  

  // APFIT muechk parameters //
  Muechk_Init();
  
  // Set size of time scanning windows
  Muechk_TimeWindowCoarse = 50; // ns
  Muechk_TimeWindowFine   = 30; // ns
  
  // Set threshold for number of hits in the scanning window
  Muechk_nThrFrac = 5; 
  Muechk_nThrMin  = 30;

  Muechk_fstBin = 5; // *10 ns
  Muechk_nBinsScanWin = 5; // *10 ns
  Muechk_prntPeakScanEnd = 40; // *10 ns
  Muechk_n50Thr = 9;

  Muechk_excessThr = 50;      // !SK2
  Muechk_signifThr = 19.9/3;  // !SK2

  // Sub-event builder parameters
  this->Set_fQuiet(TRuntimeParameters::Get().GetParameterI("spliTChan.fQuiet"));
  this->SetCluster_nBinsScanWin(TRuntimeParameters::Get().GetParameterI("spliTChan.Cluster_nBinsScanWin"));
  this->SetCluster_nHitThr(TRuntimeParameters::Get().GetParameterI("spliTChan.Cluster_nHitThr"));
  this->SetCluster_nHitThrFall(TRuntimeParameters::Get().GetParameterI("spliTChan.Cluster_nHitThrFall"));
  this->SetCluster_tBinSepThresh(TRuntimeParameters::Get().GetParameterI("spliTChan.Cluster_tBinSepThresh")); 
  //this->SetCluster_tExtendMin(TRuntimeParameters::Get().GetParameterI("spliTChan.Cluster_tExtendMin")); 

  this->SetdoTOFcorr(TRuntimeParameters::Get().GetParameterI("spliTChan.doTOFcorr")); 
  this->SetWaterRefractiveIndex(TRuntimeParameters::Get().GetParameterD("spliTChan.WaterRefractiveIndex")); 
  this->Setq_thresh(TRuntimeParameters::Get().GetParameterD("spliTChan.q_thresh"));
  this->SetMuechk_TimeWindowCoarse(TRuntimeParameters::Get().GetParameterD("spliTChan.Muechk_TimeWindowCoarse"));
  this->SetMuechk_TimeWindowFine(TRuntimeParameters::Get().GetParameterD("spliTChan.Muechk_TimeWindowFine"));
  this->SetMuechk_nThrFrac(TRuntimeParameters::Get().GetParameterI("spliTChan.Muechk_nThrFrac"));
  this->SetMuechk_nThrMin(TRuntimeParameters::Get().GetParameterI("spliTChan.Muechk_nThrMin"));
  this->SetMuechk_fstBin(TRuntimeParameters::Get().GetParameterI("spliTChan.Muechk_fstBin"));
  this->SetMuechk_nBinsScanWin(TRuntimeParameters::Get().GetParameterI("spliTChan.Muechk_nBinsScanWin"));
  this->SetMuechk_prntPeakScanEnd(TRuntimeParameters::Get().GetParameterI("spliTChan.Muechk_prntPeakScanEnd"));
  this->SetMuechk_n50Thr(TRuntimeParameters::Get().GetParameterI("spliTChan.Muechk_n50Thr"));
  this->SetMuechk_excessThr(TRuntimeParameters::Get().GetParameterD("spliTChan.Muechk_excessThr"));
  this->SetMuechk_signifThr(TRuntimeParameters::Get().GetParameterD("spliTChan.Muechk_signifThr")); 
  
  //TOF method parameters
  this->SetFixedTBeforePk(TRuntimeParameters::Get().GetParameterD("spliTChan.FixedTBeforePk"));
  this->SetFixedTAfterPk(TRuntimeParameters::Get().GetParameterD("spliTChan.FixedTAfterPk"));
  this->SetNHitThreshFact(TRuntimeParameters::Get().GetParameterD("spliTChan.nHitThreshFact"));
  if (!fQuiet) cout << "Cluster_nBinsScanWin ==== " << Cluster_nBinsScanWin  << endl;

  // Debugging histograms //
  //double t_low = -200;  // For TOF corrected h_tisk's
  //double t_high = 1500;

  double t_low = -5000;
  double t_high = 20000;

  double t_width = 10;
  int nbins = (int)((t_high-t_low)/t_width);

  if (writeHistos) {
    h_tiskz = new TH1F("h_tiskz","Hit Times (All);Time (ns)",nbins,t_low,t_high);
    
    for (int icand=0; icand<10; icand++) {
      h_tisk_corr0[icand] = new TH1F(Form("h_tisk_corr0_%d",icand),Form("Hit Times for Cluster %d;Raw Time (ns)",icand),nbins,t_low,t_high);
      h_tisk_corr1[icand] = new TH1F(Form("h_tisk_corr1_%d",icand),Form("Hit Times for Cluster %d;TOF Corrected Time (ns)",icand),nbins,t_low,t_high);
    }
    
  }
    
  t_low = -5000;
  t_high = 10000;
  t_width = 10;
  int nbinsx = (int)((t_high-t_low)/t_width);
  double q_low = -1;
  double q_high = 20;
  double q_width = 0.1;
  int nbinsy = (int)((q_high-q_low)/q_width);
  
  h_qiskz_vs_tiskz = 0;

  if (writeHistos)
    h_qiskz_vs_tiskz = new TH2F("h_qiskz_vs_tiskz","Hit Charge vs Time;Time (ns); Charge (pe)", nbinsx, t_low, t_high, nbinsy, q_low, q_high);

}

spliTChan::~spliTChan(){
  
  if (h_tiskz) h_tiskz->Delete();

  for (int icand=0; icand<MAXCAND; icand++) {
    if (h_tisk_corr0[icand]) h_tisk_corr0[icand]->Delete();
    if (h_tisk_corr1[icand]) h_tisk_corr1[icand]->Delete();
  }
  
  if (h_qiskz_vs_tiskz) h_qiskz_vs_tiskz->Delete();

  if (t_ind) delete [] t_ind;
  if (chrgz) delete [] chrgz;
  if (tHitz) delete [] tHitz;

  if (t) delete [] t;
  if (q) delete [] q;

  for (int ic=0; ic<MAXCAND; ic++) {
    if (chrg[ic]) delete [] chrg[ic];
    if (tHit[0][ic]) delete [] tHit[0][ic];
    if (tHit[1][ic]) delete [] tHit[1][ic];
  }


  std::cout << "Destroyed spliTChan" << std::endl;
}

////////////////////////////////////////////////////////
//
//  spliTChan::GetHitInfo(double *vertex)
//
//  This function copies the hit arrays(PMTall_*)
//  in fQEventManager into local arrays in this class.
//  A charge threshold of 0.3 pe is applied for usable
//  hits. The hit times are corrected for TOF given an
//  input vertex. If no vertex is passed, i.e. vertex=0,
//  then no TOF correction is applied.
//
////////////////////////////////////////////////////////
int spliTChan::GetHitInfo(double *vertex){

  
#ifndef NOSKLIBRARIES
  // Only set at first pass in case it changes
  if (!vertex) t0sk = skheadqb_.it0sk/COUNT_PER_NSEC;
#else
  if (!vertex) t0sk = 950.0/COUNT_PER_NSEC;
#endif
  
  nhitPMTz = fQEventManager::Get()->PMTall_nhit;
  
  nUse = 0;
  
  if (vertex) {
    if (!fQuiet) cout << "spliTChan::GetHitInfo(): Applying TOF correction with vertex = "
      << vertex[0] << " " << vertex[1] << " " << vertex[2] << endl;
  } else
    if (!fQuiet) cout << "spliTChan::GetHitInfo(): No TOF correction." << endl;
  
  for (int ihit=0; ihit<nhitPMTz; ihit++) {
    
    chrgz[ihit]=fQEventManager::Get()->PMTall_q[ihit];
    tHitz[ihit]=fQEventManager::Get()->PMTall_t[ihit];
    
    int icab = fQEventManager::Get()->PMTall_icab[ihit]; // C++ indexing
    
    // Ignore hits that are below some charge threshold
    if (chrgz[ihit]<=q_thresh) continue;
    
    // Ignore hits with no time info
    if (tHitz[ihit]==0) {
      if (!fQuiet) cout << "spliTChan::GetHitInfo() Warning: tiskz["
        << ihit << "] = 0" << endl;
      continue;
    }
    
    // Subtract TOF if given a vertex
    double dist=0;
    
    //if (vertex && doTOFcorr) {
    if (vertex) {
      TVector3 d;
      d[0] = fiTQun_shared::Get()->PMTpos[icab][0]-vertex[0];
      d[1] = fiTQun_shared::Get()->PMTpos[icab][1]-vertex[1];
      d[2] = fiTQun_shared::Get()->PMTpos[icab][2]-vertex[2];
      dist = d.Mag();
    }
    
    // Create array of useable hits
    //array for MakeClusters//amis
    ttc[0][nUse] = tHitz[ihit];
    ttc[1][nUse] = tHitz[ihit]-dist/cwt;
    //
    t[nUse] = tHitz[ihit]-dist/cwt;
    q[nUse] = chrgz[ihit];
    nUse++;
    
  }
  if(!nUse) return nUse;
  
  // Sort PMT hit times in ascending order
  int sortzv_mode = 1;
  int sortzv_nway = 0;
  int sortzv_nsort = 0;
  sortzv_(t, t_ind, &nUse, &sortzv_mode, &sortzv_nway, &sortzv_nsort);
  
  // Convert to C++ indexing
  for (int ihit=0; ihit<nUse; ihit++)
    t_ind[ihit] = t_ind[ihit]-1;
  
  // Initialize unused array elements
  for (int ihit=nUse; ihit<fMAXPMTZ; ihit++) t_ind[ihit] = fMAXPMTZ-1;
  t[fMAXPMTZ-1] = 1000000;
  
  if (writeHistos) GetTiskzHisto();

  return nUse;

}


void spliTChan::MuechkWrap(double *vertex, int imethod) {
  
  if (vertex) GetHitInfo(vertex);

  Muechk_Init();
    
  int i_oya = Muechk_GetParentPeakTime();

  if (i_oya < 0) {
    cout << "spliTChan::MuechkWrap() Warning: Muechk_tOya undetermined" << endl;
    return;
  }

  Muechk_FindPeakCands(i_oya);
  MuechkSubEvent();
  FillMuechkCommon(imethod);

  // Reset timing distribution to raw for writing out
  if (vertex) GetHitInfo(0);

  return;
}

void spliTChan::Muechk_Init() {

  // Initialize output variables that go to tree
  Muechk_nCand = 0;
  Muechk_tOya = 0;
  memset( Muechk_bg, 0, sizeof Muechk_bg );  
  memset( Muechk_mean, 0, sizeof Muechk_mean );  
  memset( Muechk_excess, 0, sizeof Muechk_excess );  
  memset( Muechk_signif, 0, sizeof Muechk_signif );  
  memset( Muechk_tDiff, 0, sizeof Muechk_tDiff );  
  memset( Muechk_icluster, -1, sizeof Muechk_icluster );  

}

//void spliTChan::Muechk_CommonReset() {
//
//  memset( muechk_ncand   , 0, sizeof muechk_ncand   );
//  memset( muechk_tpeak   , 0, sizeof muechk_tpeak   );
//  memset( muechk_bg      , 0, sizeof muechk_bg      );
//  memset( muechk_mean    , 0, sizeof muechk_mean    );
//  memset( muechk_excess  , 0, sizeof muechk_excess  );
//  memset( muechk_signif  , 0, sizeof muechk_signif  );
//  memset( muechk_icluster,-1, sizeof muechk_icluster);
//
//}



///////////////////////////////////////////////////////////
//
//  spliTChan::Muechk_GetParentPeakTime()
//
//  Taken from ${ATMPD_ROOT}/src/recon/muelib/muechk_gate.F
//  after line with "---------find parent peak".
//
//  Returns the starting time of the first peak of hits in 
//  a 30 ns time window.
//
///////////////////////////////////////////////////////////
int spliTChan::Muechk_GetParentPeakTime() {

  Muechk_tOya = -99999;
  
  // Set threshold for number of hits in the scanning window
  int skq_nqisk = fQEventManager::Get()->SetGate(-1000/COUNT_PER_NSEC+1000.,1496/COUNT_PER_NSEC+1000.,false,false);
  int nThr = skq_nqisk/Muechk_nThrFrac;
  if (nThr <= Muechk_nThrMin) nThr = Muechk_nThrMin;
  if (!fQuiet) cout << "spliTChan::Muechk_GetParentPeakTime(): nThr, NQISK = " << nThr << ", " << skq_nqisk << endl;

  int i=0;
  do {
  
    //cout << "fitqun " << i << " " << t[t_ind[i+nThr]] << " " << t[t_ind[i]] << " " << t[t_ind[i+nThr]]-t[t_ind[i]] << endl;
  
    i++;
    
    //if (i+nThr >= nUse) return Muechk_tOya;
    if (i+nThr >= nUse) return i;
  
  } while ( t[t_ind[i+nThr]]-t[t_ind[i]] > Muechk_TimeWindowCoarse ); 
  
  int i0 = i;
  //cout << "fitqun i0 = " << i0 << endl;
  
  int j = i0+1;
  int n_oya = 0;
  int i_oya = -1;
  do {
  
    do {
      j=j+1;
    } while ( t[t_ind[j]]-t[t_ind[i]] <= Muechk_TimeWindowFine ); // 30 ns time window
  
    int n=j-i;  // # of hits in this 30 ns time window
    if(n > n_oya) {
      n_oya = n;  // Set the most number of hits in a 30 ns time window
      i_oya = i;  // Starting index of this time window 
    }
  
    i++;
  
  } while( t[t_ind[i]]-t[t_ind[i0]] <= Muechk_TimeWindowCoarse );
  
  if (i_oya>=0) 
    Muechk_tOya = t[t_ind[i_oya]]; // <-Patrick starting time of 30 ns cluster of hits
  else 
    if (!fQuiet) cout << "spliTChan::Muechk_GetParentPeakTime() Warning: i_oya < 0" << endl;
  
  if (!fQuiet) cout << "Primary peak time (Muechk_tOya), hits (n_oya) = " << Muechk_tOya << ", " << n_oya << endl << endl;
  
  return i_oya;
}


void spliTChan::Muechk_FindPeakCands(int i_oya) {

  const int nBinMax = 3000;
  float binWidth = 10;

  TH1I h_tdiff("h_tdiff","Time Difference from Parent Peak; #Deltat (ns)", nBinMax, 0, nBinMax*binWidth);

  int use_bin[nBinMax+2] = {0};
  
  Muechk_tOya = t[t_ind[i_oya]];

  // Fill histogram with time differences from parent peak // 
  for (int ihit=0; ihit<nUse; ihit++)
    if ( t[ihit] > Muechk_tOya ) h_tdiff.Fill(t[ihit]-Muechk_tOya);
  
  //for (int ibin=1; ibin<=nBinMax; ibin++)
  //cout << "fitqun thist " << ibin << " " << h_tdiff.GetBinContent(ibin) << endl;

  // Coarse search near parent peak // 
  for (int bin=Muechk_fstBin; bin<=Muechk_prntPeakScanEnd+Muechk_nBinsScanWin; bin++) { 
    //for (int bin=Muechk_fstBin; bin<=45; bin++) { // 450 ns (original from muechk_gate.F)
    if( h_tdiff[bin] > h_tdiff[bin-1]+5 ) {
      if( h_tdiff[bin]+h_tdiff[bin+1] > h_tdiff[bin-1]+10 ) {
	float mean = (h_tdiff[bin-1]+h_tdiff[bin-2]+h_tdiff[bin-3])/3;
	float sigma = sqrt(mean);
	if( h_tdiff[bin] > mean+sigma*3 ) {
	  use_bin[bin] = 1;
	}
      }
    }
  }
 
  // Coarse search far away from parent peak // 
  int itimer = 0;
  for (int bin=Muechk_prntPeakScanEnd; bin<=nBinMax-(Muechk_nBinsScanWin-1); bin++) {
  //for (int bin=40; bin<=nBinMax-4; bin++) { // 400 ns (original from muechk_gate.F)

    int n50 = 0;
    for (int i=0; i<Muechk_nBinsScanWin; i++)
      n50 += h_tdiff[bin+i];

    int n50max, binmax, binsta;

    if (n50 > Muechk_n50Thr) {
      if (itimer == 0)  {
	itimer = 1;
	n50max = n50;
	binmax = bin;
	binsta = bin;
      } else {
	if (n50 > n50max) {
	  n50max = n50;
	  binmax = bin;
	}
      }
    } else {
      if (itimer == 1) {
	//if (bin-binsta > 4) { (original from muechk_gate.F)
	if (bin-binsta >= Muechk_nBinsScanWin) {
	  use_bin[binmax]=1;
	}
      }
      itimer = 0;
    }
    
  }

  // Do not use re-counted neighbouring bins that are probably part of the same peak
  for (int bin=Muechk_fstBin; bin<=nBinMax-1; bin++) {
    if (use_bin[bin] == 1) {
      if (use_bin[bin-2] == 1 && use_bin[bin-1] == 1 ) {
	use_bin[bin]=0;
      }
    }
  }
  for (int bin=Muechk_fstBin; bin<=nBinMax-1; bin++) {
    if (use_bin[bin] == 1) {
      if (use_bin[bin-1] == 1 ) {
	use_bin[bin]=0;
      }
    }
  }
  
  //cout << "spliTChan::Muechk_FindPeakCands(): Candidate bins from coarse search" << endl;
  //cout << "bin  bin_low_edge" << endl;
  //for (int bin=1; bin<=nBinMax; bin++)
  //  if (use_bin[bin] == 1) 
  //    cout << bin << " " << h_tdiff.GetBinLowEdge(bin) << endl;
  //cout << endl;

  // Fine search //

  //cout << "spliTChan::Muechk_FindPeakCands(): Starting fine search..." << endl;

  int i0 = i_oya;
  for (int bin=1; bin<=Muechk_fstBin-1; bin++)
    i0 = i0+h_tdiff[bin];
    
  Muechk_nCand = 0;
  for (int bin=Muechk_fstBin; bin<=nBinMax-1; bin++) {
    
    if (use_bin[bin] == 1) {

      // find peak and count nhit in window //
      int n_peak = 0;
      int i_peak = -1;
      float t_peak;

      int i0_min = i0 - h_tdiff[bin-1];
      int i0_max = i0 + h_tdiff[bin] + h_tdiff[bin+1];

      //cout << "fitqun blah " << bin << " " << i0 << " " << h_tdiff.GetBinContent(bin-1) << " " << h_tdiff.GetBinContent(bin) << " " << h_tdiff.GetBinContent(bin+1) << endl;
      //cout << "fitqun blah " << bin << " " << i0_min << " " << i0_max << " " << nUse << endl;
    
      if (i0_max >= nUse-1) i0_max = nUse-2;

      for (int i=i0_min; i<=i0_max; i++) {
      
	int n=1;
	do { 
	  n++;
	} while (t[t_ind[i+n]] - t[t_ind[i]] <= Muechk_TimeWindowFine );
		
	if (n > n_peak) {
	  n_peak = n;
	  i_peak = i;
	  t_peak = t[t_ind[i_peak]];
	}
      }

      Muechk_tDiff[Muechk_nCand] = t_peak - Muechk_tOya;
      //cout << "bin = " << bin << endl;
      //cout << "t_peak, n_peak, tDiff = " << t_peak << ", " << n_peak << ", " << Muechk_tDiff[Muechk_nCand] << endl;

      // estimate background //
      double tWinBg;

      // Maybe this can be tweaked
      if (Muechk_tDiff[Muechk_nCand] <= 70) tWinBg = 20;
      else if (Muechk_tDiff[Muechk_nCand] <= 300) tWinBg = 30;
      else if (Muechk_tDiff[Muechk_nCand] <= 350) tWinBg = 50;
      else tWinBg = 100;

      if (Muechk_nCand > 0) {
	if (Muechk_tDiff[Muechk_nCand] - Muechk_tDiff[Muechk_nCand-1] <= tWinBg + Muechk_TimeWindowFine) {
	  tWinBg = Muechk_tDiff[Muechk_nCand] - Muechk_tDiff[Muechk_nCand-1] - Muechk_TimeWindowFine - 10; // Maybe this can be tweaked
	  if (tWinBg < 20) tWinBg = 20; // Maybe this can be tweaked
	}
      }

      Muechk_bg[Muechk_nCand] = 0;

      ///****************************
      ///*****************************
      /// CRASH IN BELOW WHILE LOOP (in t_ind.at() overflow?)
      do {
	Muechk_bg[Muechk_nCand]++;
	//std::cout<<"Muechk nCand="<<Muechk_nCand
	//	 <<" i_peak="<<i_peak
	//		 <<" bg="<<Muechk_bg[Muechk_nCand]
	//		 <<" i_p-bg-ncand="<<i_peak-Muechk_bg[Muechk_nCand]-1
	  //<<" t_ind.size="<<t_ind.size()
	//	 <<std::endl;
	// fail here:
	if ( i_peak-Muechk_bg[Muechk_nCand]-1 < 0 || i_peak <0 ) return;
      } while( t[t_ind[i_peak]] - t[t_ind[i_peak-Muechk_bg[Muechk_nCand]-1]] <= tWinBg );
      Muechk_mean[Muechk_nCand] = Muechk_bg[Muechk_nCand] / tWinBg * Muechk_TimeWindowFine;
      float sigma = sqrt(Muechk_mean[Muechk_nCand]);
      if (sigma <= 3) sigma = 3;

      //cout << "tWinBg, n_bg, mean, sigma = " << tWinBg << ", " << Muechk_bg[Muechk_nCand] << ", " << Muechk_mean[Muechk_nCand] << ", " << sigma << endl;

      // compare S/N //
      Muechk_excess[Muechk_nCand] = n_peak - Muechk_mean[Muechk_nCand];
      Muechk_signif[Muechk_nCand] = Muechk_excess[Muechk_nCand] / sigma;
      
      int iflag = 0;
      if (Muechk_nCand > 0) {
	if (Muechk_tDiff[Muechk_nCand]-Muechk_tDiff[Muechk_nCand-1] < 1000) { // Maybe this can be tweaked
	  if (Muechk_excess[Muechk_nCand] < 10) iflag = -1; // Maybe this can be tweaked
	  iflag = 1;
	}
      }
      
      if (iflag != -1) {
	if ( (Muechk_tDiff[Muechk_nCand] > 600 && iflag == 0) // Maybe this can be tweaked
	     || Muechk_excess[Muechk_nCand] >= Muechk_excessThr
	     || Muechk_signif[Muechk_nCand] > Muechk_signifThr) {
	  
	  Muechk_nCand++;
	
	  // Find candidates that are in the same 30 ns time window	
	  if (Muechk_nCand > 1) {
	    if (Muechk_tDiff[Muechk_nCand-1]-Muechk_tDiff[Muechk_nCand-2] <= Muechk_TimeWindowFine) {
	    
	      Muechk_nCand--; // Probably the same decay-e

	      // Set candidate to the lower of the two significances
	      if (Muechk_signif[Muechk_nCand-1] < Muechk_signif[Muechk_nCand]) {
		Muechk_tDiff[Muechk_nCand-1]  = Muechk_tDiff[Muechk_nCand];
		Muechk_excess[Muechk_nCand-1] = Muechk_excess[Muechk_nCand];
		Muechk_signif[Muechk_nCand-1] = Muechk_signif[Muechk_nCand];
	      }
	    }
	  }
	
	  if (Muechk_nCand >= MAXCAND-1) break;
	}
      }
      // 30 continue (in muechk_gate.F)
    }
    i0 = i0 + h_tdiff[bin];
  }
  if (!fQuiet) cout << endl;

  // Convert back to absolute time of TISKZ
  for (int i=0; i<Muechk_nCand; i++) 
    Muechk_tDiff[i] += Muechk_tOya;
  
  // Shift all peak candidates in array by 1
  if (Muechk_nCand>0) {
    for (int i=Muechk_nCand-1; i>=0; i--) {
      Muechk_bg[i+1]     = Muechk_bg[i]    ;    
      Muechk_mean[i+1]   = Muechk_mean[i]  ;  
      Muechk_tDiff[i+1]  = Muechk_tDiff[i] ; 
      Muechk_excess[i+1] = Muechk_excess[i];
      Muechk_signif[i+1] = Muechk_signif[i];
    }
  }
  
  // Add parent peak to beginning of array
  int i=0;
  Muechk_bg[i]     = -1;    
  Muechk_mean[i]   = -1;  
  Muechk_tDiff[i]  = Muechk_tOya;
  Muechk_excess[i] = -1;
  Muechk_signif[i] = -1;  
  Muechk_nCand++;
  
  if (!fQuiet) {
    cout << "spliTChan::Muechk_FindPeakCands(): Fine search results" << endl;
    cout << "Muechk_nCand = " << Muechk_nCand << endl;
    for (int i=0; i<Muechk_nCand; i++) {
      cout << "i, dt (ns), excess, signif = " << i << ", " << Muechk_tDiff[i] << ", " << Muechk_excess[i] << ", " << Muechk_signif[i] << endl;
    }
    cout << endl;
  }
  

}

// Copied from muelib/muechk.F, after call to muechk_gate
void spliTChan::MuechkSubEvent() {

  float ptime[MAXCAND] = {0};

  //Count sub-events as peaks that occur 600 ns after primary peak
  int n_sub = 0;
  for (int imue=1; imue<Muechk_nCand; imue++) { // start at 1 to ignore parent peak
    if (Muechk_tDiff[imue]-Muechk_tOya>600) {
      ptime[n_sub] = Muechk_tDiff[imue];
      n_sub++;
    }
  }

  //cout << "spliTChan: n_sub = " << n_sub << endl;

  if (!n_sub) return;
  
  int Muechk_nCand_inGate = Muechk_nCand - n_sub; 
  // Muechk_nCand is now the number of "in-gate" decay-e

  for (int isub=0; isub<n_sub; isub++) {

    // First sub-event and at least 1 in-gate peak &&
    // Time difference between first sub-event and last in-gate peak 
    //if (isub==0 && Muechk_nCand_inGate>1)  // >1 to ignore parent peak
    //  cout << "spliTChan: ptime[0] diffs = " << isub << " " << ptime[0]-Muechk_tDiff[Muechk_nCand_inGate-1] << endl;
    //
    //if (isub<n_sub-1)
    //  cout << "spliTChan: ptime[isub+1] diffs = " << isub << " " << ptime[isub+1]-ptime[isub] << endl;
    //  
    //if (isub>0)
    //  cout << "spliTChan: ptime[isub-1] diffs = " << isub << " " << ptime[isub]-ptime[isub-1] << endl;

    if ( (isub==0 && Muechk_nCand_inGate>1 && ptime[0]-Muechk_tDiff[Muechk_nCand_inGate-1]<200) ||
	 //Time difference between this sub-event and next
	 (isub<n_sub-1 && ptime[isub+1]-ptime[isub]<500) ||

	 // Time difference between this sub-event and previous
	 (isub>0 && ptime[isub]-ptime[isub-1]<500) ) {

      //cout << "spliTChan: isub, ptime, excess = " << isub << " " <<  ptime[isub]-Muechk_tOya << " " << Muechk_excess[isub+1] << endl;
      
      if(isub<n_sub-1) {

	// Time difference between this sub-event and next; and hits above BG
	// (If sub-events are too close together with insufficient hits, do not 
	// count as additional Muechk_nCand decay-e)
	//if (ptime[isub+1]-ptime[isub]<500 && Muechk_excess[isub+Muechk_nCand_inGate]<25) RemoveMuechk(isub+Muechk_nCand_inGate);
	if (ptime[isub+1]-ptime[isub]<500 && Muechk_excess[isub+1]<25) RemoveMuechk(isub+Muechk_nCand_inGate); // Is this a bug? I think it should be as in the previous line
	// Anyways keep the same as original muechk for consistency (also +1 to ignore parent peak)
      }
    }
  }
  
}

void spliTChan::RemoveMuechk(int amue) {
  
  if (!fQuiet) cout << "spliTChan::RemoveMuechk(): Removing peak " << amue << " " << Muechk_excess[amue] << endl;

  for (int imue=amue; imue<Muechk_nCand-1; imue++) {
    Muechk_bg[imue] = Muechk_bg[imue+1];
    Muechk_mean[imue] = Muechk_mean[imue+1];
    Muechk_excess[imue] = Muechk_excess[imue+1];
    Muechk_signif[imue] = Muechk_signif[imue+1];
    Muechk_tDiff[imue] = Muechk_tDiff[imue+1];
  }

  Muechk_bg[Muechk_nCand-1]     = 0;
  Muechk_mean[Muechk_nCand-1]   = 0;
  Muechk_excess[Muechk_nCand-1] = 0;
  Muechk_signif[Muechk_nCand-1] = 0;
  Muechk_tDiff[Muechk_nCand-1]  = 0;
  Muechk_nCand--;
}

// Copied from Muechk_FindPeakCands but removed everything about the parent peak
void spliTChan::Muechk_FindPeakCandsFine(int i_oya) {

  const int nBinMax = 3000;
  float binWidth = 10;

  TH1I h_tdiff("h_tdiff","Time Difference from Parent Peak; #Deltat (ns)", nBinMax, 0, nBinMax*binWidth);

  int use_bin[nBinMax+2] = {0};

  Muechk_tOya = t[t_ind[i_oya]];

  // Fill histogram with time differences from parent peak // 
  for (int ihit=0; ihit<nUse; ihit++)
    if ( t[ihit] > Muechk_tOya ) h_tdiff.Fill(t[ihit]-Muechk_tOya);
  
  // Coarse search far away from parent peak // 
  int itimer = 0;
  for (int bin=1; bin<=nBinMax-(Muechk_nBinsScanWin-1); bin++) {

    int n50 = 0;
    for (int i=0; i<Muechk_nBinsScanWin; i++)
      n50 += h_tdiff[bin+i];

    int n50max, binmax, binsta;

    if (n50 > Muechk_n50Thr) {
      if (itimer == 0)  {
	itimer = 1;
	n50max = n50;
	binmax = bin;
	binsta = bin;
      } else {
	if (n50 > n50max) {
	  n50max = n50;
	  binmax = bin;
	}
      }
    } else {
      if (itimer == 1) {
	//if (bin-binsta > 4) { (original from muechk_gate.F)
	if (bin-binsta >= Muechk_nBinsScanWin) {
	  use_bin[binmax]=1;
	}
      }
      itimer = 0;
    }
    
  }

  // Do not use re-counted neighbouring bins that are probably part of the same peak
  for (int bin=1; bin<=nBinMax-1; bin++) {
    if (use_bin[bin] == 1) {
      if (use_bin[bin-2] == 1 && use_bin[bin-1] == 1 ) {
	use_bin[bin]=0;
      }
    }
  }
  for (int bin=1; bin<=nBinMax-1; bin++) {
    if (use_bin[bin] == 1) {
      if (use_bin[bin-1] == 1 ) {
	use_bin[bin]=0;
      }
    }
  }
  
  //cout << "spliTChan::Muechk_FindPeakCandsFine(): Candidate bins from coarse search" << endl;
  //cout << "bin  bin_low_edge" << endl;
  //for (int bin=1; bin<=nBinMax; bin++)
  //  if (use_bin[bin] == 1) 
  //    cout << bin << " " << h_tdiff.GetBinLowEdge(bin) << endl;
  //cout << endl;

  // Fine search //

  //cout << "spliTChan::Muechk_FindPeakCandsFine(): Starting fine search..." << endl;

  int i0 = i_oya;
  
  Muechk_nCand = 0;
  for (int bin=1; bin<=nBinMax-1; bin++) {
    
    if (use_bin[bin] == 1) {

      // find peak and count nhit in window //
      int n_peak = 0;
      int i_peak = -1;
      float t_peak;

      int i0_min = i0 - h_tdiff[bin-1];
      int i0_max = i0 + h_tdiff[bin] + h_tdiff[bin+1];
    
      if (i0_max >= nUse-1) i0_max = nUse-2;

      for (int i=i0_min; i<=i0_max; i++) {
      
	int n=1;
	do { 
	  n++;
	} while (t[t_ind[i+n]] - t[t_ind[i]] <= Muechk_TimeWindowFine );
		
	if (n > n_peak) {
	  n_peak = n;
	  i_peak = i;
	  t_peak = t[t_ind[i_peak]];
	}
      }
      
      Muechk_tDiff[Muechk_nCand] = t_peak; //- Muechk_tOya;
      
      Muechk_nCand++;      
	
      if (Muechk_nCand >= MAXCAND) break;
      
    }
    i0 = i0 + h_tdiff[bin];
  }
  
  // Convert back to absolute time of TISKZ
  //for (int i=0; i<Muechk_nCand; i++) 
  //Muechk_tDiff[i] += Muechk_tOya;

  // (something weird here where sometimes parent peak is found so don't refill)
  bool parentFound = 0;
  for (int ipeak=0; ipeak<Muechk_nCand; ipeak++)
    if (fabs(Muechk_tDiff[ipeak]-Muechk_tOya) < 0.00001)
      parentFound = 1;
  
  // Shift all peak candidates in array by 1
  if (!parentFound) {
        
    for (int ipeak=Muechk_nCand-1; ipeak>=0; ipeak--) {
      Muechk_bg[ipeak+1]     = Muechk_bg[ipeak]    ;    
      Muechk_mean[ipeak+1]   = Muechk_mean[ipeak]  ;  
      Muechk_tDiff[ipeak+1]  = Muechk_tDiff[ipeak] ; 
      Muechk_excess[ipeak+1] = Muechk_excess[ipeak];
      Muechk_signif[ipeak+1] = Muechk_signif[ipeak];
    }

    // Add parent peak to beginning of array
    int ipeak=0;
    Muechk_bg[ipeak]     = -1;    
    Muechk_mean[ipeak]   = -1;  
    Muechk_tDiff[ipeak]  = Muechk_tOya;
    Muechk_excess[ipeak] = -1;
    Muechk_signif[ipeak] = -1;  
    Muechk_nCand++;

  }

  if (!fQuiet) {
    cout << "spliTChan::Muechk_FindPeakCandsFine(): Fine search results" << endl;
    cout << "Muechk_nCand = " << Muechk_nCand << endl;
    for (int i=0; i<Muechk_nCand; i++) {
      cout << "i, dt (ns), excess, signif = " << i << ", " << Muechk_tDiff[i] << ", " << Muechk_excess[i] << ", " << Muechk_signif[i] << endl;
    }
    cout << endl;
  }
  

}


void spliTChan::Cluster_Init() {

  // Initialize output variables that go to tree
  Cluster_nCand = 0;
  memset( Cluster_tStart, 0, sizeof Cluster_tStart );  
  memset( Cluster_tEnd, 0, sizeof Cluster_tEnd );  
  memset( Cluster_vertex, 0, sizeof Cluster_vertex );  
  memset( Cluster_totQ, 0, sizeof Cluster_totQ );  
  memset( Cluster_nHits, 0, sizeof Cluster_nHits );  
  //memset( Cluster_goodflag, 1, sizeof Cluster_goodflag );  
  //memset( Cluster_nPeaks, 0, sizeof Cluster_nPeaks );
//  memset( chrg, 0, sizeof chrg );    
//  memset( tHit, 0, sizeof tHit );

  for (int icand=0; icand<MAXCAND; icand++) {

  memset( chrg[icand], 0, sizeof chrg[icand] );    
  memset( tHit[0][icand], 0, sizeof tHit[0][icand] );
  memset( tHit[1][icand], 0, sizeof tHit[1][icand] );
    Cluster_goodflag[icand] = 1;
    if (h_tisk_corr0[icand]) h_tisk_corr0[icand]->Reset();
    if (h_tisk_corr1[icand]) h_tisk_corr1[icand]->Reset();
  }
  
  Cluster_timeOfPeak.clear();
  Cluster_ipeak.clear();

}


void spliTChan::Cluster_Search(double *vertex) {

  GetHitInfo(vertex);
    
  Cluster_Init();

  const int nBinMax = 100000;
  float binWidth = 10;

  TH1I h_tdiff("h_tdiff","Time Difference from Parent Peak; #Deltat (ns)", nBinMax, 0, nBinMax*binWidth);

  int use_bin[nBinMax+2] = {0};

  int i_oya = 0;
  float Cluster_tOya = t[t_ind[i_oya]];

  // Fill histogram with time differences from parent peak // 
  for (int ihit=0; ihit<nUse; ihit++)
    if ( t[ihit] > Cluster_tOya ) h_tdiff.Fill(t[ihit]-Cluster_tOya);
  
  //for (int ibin=1; ibin<=nBinMax; ibin++)
  //cout << "fitqun thist " << ibin << " " << h_tdiff.GetBinContent(ibin) << endl;

  if (h_tdiff[0]) {
    cout << "spliTChan::Cluster_Search() Error: h_tdiff underflow = " << h_tdiff[0] << endl;
    exit (-1);
  }

  if (h_tdiff[nBinMax+1]) {
    cout << "spliTChan::Cluster_Search() Error: h_tdiff overflow = " << h_tdiff[nBinMax+1] << endl;
    exit (-1);
  }

  int itimer = 0;
  vector<int> use_binsta;
  vector<int> use_binend;
  Cluster_nCand = 0;
  
  int binend;
  for (int bin=1; bin<=nBinMax-(Cluster_nBinsScanWin-1); bin++) {

    int nHit = (int)h_tdiff.Integral(bin,bin+Cluster_nBinsScanWin-1);
    
    if (itimer==0 && nHit > Cluster_nHitThr) {
      itimer = 1;
      use_binsta.push_back(bin);

      // Use bin at the end of scanning window (+1 for high edge)
      binend = bin+Cluster_nBinsScanWin+1;
    }
    else if (itimer==1 && nHit > Cluster_nHitThrFall) {

      // Use bin at the end of scanning window (+1 for high edge)
      binend = bin+Cluster_nBinsScanWin+1;

    } else {
      if (itimer == 1) {
	if (bin-use_binsta.back() >= Cluster_nBinsScanWin) {
	  
	  if (!fQuiet) cout << "SubEvent " << Cluster_nCand << " startBin, startTime  = " << use_binsta.back() << ", " << h_tdiff.GetBinLowEdge(use_binsta.back())+Cluster_tOya << endl;

	  // Extend end of time window to catch stray (reflected) hits 
	  // that may lie under noise (linear based on plot of furthest missed time
	  // vs # of hits in the cluster, and each event weighted by % of missed hits)
	  int nHitsInCluster = (int)h_tdiff.Integral(use_binsta.back(), binend-1);
	  //float timeExtension = TMath::Max(0., (float)(-nHitsInCluster*300./500. + 300)); // 1st iteration
	  //float timeExtension = TMath::Max(Cluster_tExtendMin, (float)(-nHitsInCluster*320./2000. + 320));
	  float timeExtension = TMath::Max(0., (-nHitsInCluster*320./800. + 320));

	  binend += (int)ceil(timeExtension/h_tdiff.GetBinWidth(binend));

	  // Extend low edge by 30 ns to make sure we collect all hits
	  // (also FQ peak finder seems to be too early if low charge)
	  use_binsta.at(use_binsta.size()-1) -= 3;

	  if (timeExtension)
	    if (!fQuiet) cout << "Extending time window by " << timeExtension << "(" << binend << ") for " << nHitsInCluster << " hit cluster." << endl;

	  // Ignore clusters that start too close to previous one
	  if ( (Cluster_nCand > 0) && 
	       (use_binsta.back()-use_binend[Cluster_nCand-1] <= Cluster_tBinSepThresh) ) {
	    
	    // Remove current cluster start time
	    if (!fQuiet) cout << "Ignoring this SubEvent since it is too close to previous one." << endl;
	    use_binsta.pop_back();
	    
	    // Remove previous cluster start time
	    if (!fQuiet) cout << "Extending previous SubEvent window to include this cluster." << endl;
	    use_binend.pop_back();

	    if (!fQuiet) cout << "SubEvent " << Cluster_nCand-1 << " endBin, endTime = " << binend  << ", " << h_tdiff.GetBinLowEdge(binend)+Cluster_tOya << endl;
	    use_binend.push_back(binend);
	    
	  }
	  
	  else { // This is a good cluster
	    if (!fQuiet) cout << "SubEvent " << Cluster_nCand << " endBin, endTime = " << binend  << ", " << h_tdiff.GetBinLowEdge(binend)+Cluster_tOya << endl;
	    use_binend.push_back(binend);
	    Cluster_nCand++;
	  }

	} else {
	  use_binsta.pop_back();
	}

      }
      itimer = 0;
    }
  }
  
  if (use_binsta.size() != Cluster_nCand) {
    cout << "binsta size != nCand " << use_binsta.size() << " != " <<  Cluster_nCand << endl;
    exit (-1);
  }
  
  if (use_binend.size() != Cluster_nCand) {
    cout << "binend size != nCand " << use_binend.size() << " != " <<  Cluster_nCand << endl;
    exit (-1);
  }

  for (int icluster=0; icluster<Cluster_nCand; icluster++) {
    Cluster_tStart[icluster] = h_tdiff.GetBinLowEdge(use_binsta[icluster])+Cluster_tOya;
    Cluster_tEnd[icluster] = h_tdiff.GetBinLowEdge(use_binend[icluster])+Cluster_tOya;
  }

  RemoveOutGateClusters();
  CalcClusterProperties();
  //RemoveLowHitFirstClusters(); // Apply peak finding rejection instead (currently in runfiTQun.cc)

  FillClusterCommon();
  
  if (writeHistos) {
    FillTQbanksAll();
  }
  
}

// Remove too far outside original in-gate time window (typically -5 us, 35us)
void spliTChan::RemoveOutGateClusters() {
  
  for (int icluster=0; icluster < Cluster_nCand; icluster++) {
    
    int skq_nqisk = fQEventManager::Get()->SetGate(this->Cluster_tStart[icluster],this->Cluster_tEnd[icluster],true,false);
    
    if (skq_nqisk<=0) {
      if (!fQuiet) cout << "runfiTQun Warning: skq_nqisk = " << skq_nqisk << " for icluster = " << icluster << ". Removing cluster [" << this->Cluster_tStart[icluster] << ", " << this->Cluster_tEnd[icluster] << "]" << endl;
      this->RemoveCluster(icluster); // This shifts Cluster_* array variables down by one and decrements Cluster_nCand
      icluster--; 
      
    }
  }
}


void spliTChan::CalcClusterProperties() {
  
  for (int icluster=0; icluster < Cluster_nCand; icluster++) {
    Cluster_nHits[icluster] = 0;
    Cluster_totQ[icluster]  = 0;
    for (int ihit=0; ihit<nUse; ihit++) {

      if (Cluster_tStart[icluster]<t[ihit] && t[ihit]<Cluster_tEnd[icluster]) {
        Cluster_nHits[icluster]++;
        Cluster_totQ[icluster] += q[ihit];
      }
    }
  }
  
}


// Attempt to cut out all false first clusters since fiTQun currently 
// only runs 1R and pi0 stuff on first cluster
void spliTChan::RemoveLowHitFirstClusters() {

  // This should be a global parameter, but seg. faults
  // when added to spliTChan.h
  //float Cluster_totQFrstClstrCut = 0.25;  // per ns (normalized to time window), estimated from atmospheric data
  float Cluster_nHitsFrstClstrCut = 0.115;  // per ns (normalized to time window), estimated from atmospheric data

  for (int icluster=0; icluster<Cluster_nCand; icluster++) {
  
    float hitsovert = Cluster_nHits[icluster]/(Cluster_tEnd[icluster]-Cluster_tStart[icluster]);

    //if (Cluster_totQ[icluster]/(Cluster_tEnd[icluster]-Cluster_tStart[icluster]) < Cluster_totQFrstClstrCut) {
    if (hitsovert < Cluster_nHitsFrstClstrCut) {
    if (!fQuiet) cout << "spliTChan::RemoveLowHitFirstClusters(): Tagging bad cluster " << icluster << "nHits/T = " << hitsovert << " < " << Cluster_nHitsFrstClstrCut << endl;
      Cluster_goodflag[icluster] = 0;
    }
    else break;
    
  }
  
}

// Tag fake clusters as those without a coresponding peak
void spliTChan::RemoveNoPeakClusters(int imethod) {
  
  for (int icluster=0; icluster<Cluster_nCand; icluster++) {
    if (!fitqunse_.cluster_npeaks[icluster][imethod]) {
      if (!fQuiet) cout << "spliTChan::RemoveNoPeakClusters(): Tagging bad cluster " << icluster << " with no peaks." << endl;
      Cluster_goodflag[icluster] = 0;
    }
  }

  FillClusterCommon();

}


void spliTChan::Cluster_AssociateMuechk() {
  
  Cluster_timeOfPeak.clear();
  Cluster_ipeak.clear();
  for (int icluster=0; icluster<Cluster_nCand; icluster++) {
      
    std::vector<std::vector<double> > timeOfPeaksTempTemp;
    std::vector<std::vector<int> > ipeakTempTemp;

    for (int imethod=0; imethod<NPEAKFINDS; imethod++) {

      // Count peaks in cluster 
      fitqunse_.cluster_npeaks[icluster][imethod]=0;
      std::vector<double> timeOfPeaksTemp;
      std::vector<int> ipeakTemp;      
      for (int ipeak=0; ipeak<fitqunse_.muechk_ncand[imethod]; ipeak++) {
	
	if (Cluster_tStart[icluster]<=fitqunse_.muechk_tpeak[imethod][ipeak] &&
	    fitqunse_.muechk_tpeak[imethod][ipeak]<=Cluster_tEnd[icluster]) {
	
	  fitqunse_.muechk_icluster[imethod][ipeak] = icluster;

	  ipeakTemp.push_back(ipeak);
	  timeOfPeaksTemp.push_back(fitqunse_.muechk_tpeak[imethod][ipeak]);
	  fitqunse_.cluster_npeaks[icluster][imethod]++;
	}
      }
      ipeakTempTemp.push_back(ipeakTemp);
      timeOfPeaksTempTemp.push_back(timeOfPeaksTemp); 

      // Count peaks in cluster 
      //Cluster_nPeaks[icluster]=0;
      //std::vector<double> timeOfPeaksTemp;
      //std::vector<double> ipeakTemp;      
      //for (int ipeak=0; ipeak<Muechk_nCand; ipeak++) {
      //	
      //	if (Cluster_tStart[icluster]<=Muechk_tDiff[ipeak] &&
      //	    Muechk_tDiff[ipeak]<=Cluster_tEnd[icluster]) {
      //	
      //	  Muechk_icluster[ipeak] = icluster;
      //
      //	  ipeakTemp.push_back(ipeak);
      //	  timeOfPeaksTemp.push_back(Muechk_tDiff[ipeak]);
      //	  Cluster_nPeaks[icluster]++;
      //	}
      //}
      //Cluster_ipeak.push_back(ipeakTemp);
      //Cluster_timeOfPeak.push_back(timeOfPeaksTemp); 

    }
    
    Cluster_ipeak.push_back(ipeakTempTemp);
    Cluster_timeOfPeak.push_back(timeOfPeaksTempTemp); 

  }

  FillClusterPeakCommon();
}


void spliTChan::FillClusterCommon() {

  if (Cluster_nCand>MAXSE) {
    cout << "spliTChan::FillClusterCommon() Warning: Number of clusters exceeded " << MAXSE << endl;
   Cluster_nCand = MAXSE;    
  }
  fitqunse_.cluster_ncand = Cluster_nCand;
  
  for (int icluster=0; icluster<Cluster_nCand; icluster++) {
    fitqunse_.cluster_tstart[icluster] = Cluster_tStart[icluster];
    fitqunse_.cluster_tend[icluster] = Cluster_tEnd[icluster];
    fitqunse_.cluster_nhits[icluster] = Cluster_nHits[icluster];
    fitqunse_.cluster_totq[icluster] = Cluster_totQ[icluster];
    fitqunse_.cluster_goodflag[icluster] = Cluster_goodflag[icluster];
  }
}

void spliTChan::CopyPeakInfo(int imethod, int npeaks, float *peakt0, float *peakiness) {

  Muechk_Init();

  Muechk_nCand = npeaks;
  for (int ipeak=0; ipeak<npeaks; ipeak++) {
    Muechk_tDiff[ipeak] = peakt0[ipeak];
    Muechk_signif[ipeak] = peakiness[ipeak];
  }
  
  FillMuechkCommon(imethod);

}

void spliTChan::FillMuechkCommon(int imethod) {
  if (Muechk_nCand>MAXSE) {
    cout << "spliTChan::FillMuechkCommon() Warning: Number of peaks exceeded " << MAXSE << endl;
    Muechk_nCand = MAXSE;    
  }
  
  fitqunse_.muechk_ncand[imethod] = Muechk_nCand;
  fitqunse_.muechk_toya[imethod] = Muechk_tOya;
  
  for (int ipeak=0; ipeak<Muechk_nCand; ipeak++) {
    fitqunse_.muechk_tpeak[imethod][ipeak] = Muechk_tDiff[ipeak];
    fitqunse_.muechk_bg[imethod][ipeak] = Muechk_bg[ipeak];
    fitqunse_.muechk_mean[imethod][ipeak] = Muechk_mean[ipeak];
    fitqunse_.muechk_excess[imethod][ipeak] = Muechk_excess[ipeak];
    fitqunse_.muechk_signif[imethod][ipeak] = Muechk_signif[ipeak];
  }
}

void spliTChan::FillClusterPeakCommon() {

  //for (int ipeak=0; ipeak<Muechk_nCand; ipeak++) {
  //  fitqunse_.muechk_icluster[ipeak] = Muechk_icluster[ipeak];
  //}

  for (int icluster=0; icluster<Cluster_nCand; icluster++) {

    for (int imethod=0; imethod<NPEAKFINDS; imethod++) {

      //fitqunse_.cluster_npeaks[icluster][imethod] = Cluster_nPeaks[icluster][imethod];
    
      for (int ipeak=0; ipeak<fitqunse_.cluster_npeaks[icluster][imethod]; ipeak++) {
	fitqunse_.cluster_ipeak[icluster][imethod][ipeak] = Cluster_ipeak[icluster][imethod][ipeak];
	fitqunse_.cluster_timeofpeak[icluster][imethod][ipeak] = Cluster_timeOfPeak[icluster][imethod][ipeak];
      }
    }
  }
}

int spliTChan::FillTQbanksAll() {
  
  // It's all relative to first hit cluster which is defined immediately after
  // software trigger, but before Muechk (which then re-finds it...)
  for (int icluster=0; icluster<Cluster_nCand; icluster++) {
    
    fQEventManager::Get()->SetGate(Cluster_tStart[icluster],Cluster_tEnd[icluster],true,false);
    
    float *vertex=0;
    if (Cluster_vertex[icluster][0] && Cluster_vertex[icluster][1] && Cluster_vertex[icluster][2]) {
      
      vertex = Cluster_vertex[icluster];
      if (!fQuiet) cout << "spliTChan::FillTQbanksAll(): Applying TOF correction with vertex = "
        << vertex[0] << " " << vertex[1] << " " << vertex[2] << endl;
    }
    //else
    //cout << "spliTChan::FillTQbanksAll(): No TOF correction." << endl;
    
    for (int icab=0; icab<fNPMT; icab++) {
      
      if (fQEventManager::Get()->PMT_iHit[icab]==-1) continue;
      
      if (fQEventManager::Get()->PMT_q[icab]<=q_thresh) continue;
      
      double dist=0;
      if (vertex) {
        TVector3 d;
        d[0] = fiTQun_shared::Get()->PMTpos[icab][0]-vertex[0];
        d[1] = fiTQun_shared::Get()->PMTpos[icab][1]-vertex[1];
        d[2] = fiTQun_shared::Get()->PMTpos[icab][2]-vertex[2];
        dist = d.Mag();
      }
      
      chrg[icluster][icab] = fQEventManager::Get()->PMT_q[icab];
      tHit[(vertex!=0)][icluster][icab] = fQEventManager::Get()->PMT_t[icab]-dist/cwt;
      
      if (!vertex) h_tisk_corr0[icluster]->Fill(fQEventManager::Get()->PMT_t[icab]);
      else h_tisk_corr1[icluster]->Fill(fQEventManager::Get()->PMT_t[icab]-dist/cwt);
      
    }
  }
  
  return 0;
  
}

TH1F* spliTChan::GetTiskzHisto() {

  if (!writeHistos) return 0;

  h_tiskz->Reset();
  
  for (int ihit=0; ihit<nUse; ihit++) {
    h_tiskz->Fill(t[ihit]);
  }
  
  return h_tiskz;
  
}


void spliTChan::SetClusterVertex(int icluster, double *vtxtmp) {
  
  for (int i=0; i<nClusterVertexVars; i++) {
    Cluster_vertex[icluster][i] = vtxtmp[i];
  }
   
}


void spliTChan::RemoveCluster(int acluster) {
  
  for (int icluster=acluster; icluster<Cluster_nCand-1; icluster++) {
    Cluster_tStart[icluster] = Cluster_tStart[icluster+1];
    Cluster_tEnd[icluster] = Cluster_tEnd[icluster+1];
    Cluster_totQ[icluster] = Cluster_totQ[icluster+1];
    Cluster_nHits[icluster] = Cluster_nHits[icluster+1];
    Cluster_goodflag[icluster] = Cluster_goodflag[icluster+1];
    
    for (int i=0; i<nClusterVertexVars; i++) Cluster_vertex[icluster][i] = Cluster_vertex[icluster+1][i];
  }
  Cluster_tStart[Cluster_nCand-1] = 0;
  Cluster_tEnd[Cluster_nCand-1] = 0;
  Cluster_totQ[Cluster_nCand-1] = 0;
  Cluster_nHits[Cluster_nCand-1] = 0;
  Cluster_goodflag[Cluster_nCand-1] = 1;
  for (int i=0; i<nClusterVertexVars; i++) Cluster_vertex[Cluster_nCand-1][i] = 0;
  Cluster_nCand--;
}


TH2F* spliTChan::GetQiskzVsTiskzHisto() {

  if (!writeHistos) return 0;

  h_qiskz_vs_tiskz->Reset();
  
  for (int ihit=0; ihit<nhitPMTz; ihit++) {
    h_qiskz_vs_tiskz->Fill( tHitz[ihit], chrgz[ihit]);
  }
  
  return h_qiskz_vs_tiskz;
  
}

//new from amis//


  
int spliTChan::isValidHit(float thit, float trange){
  //make sure there are cluster hits in nint 20ns intervals of a hit
  //loop over all hits in cluster and make sure
//  cout<<"called is Valid Hit"<<endl; 
  int nhits=0;
  
 // float nhit_thresh = TMath::Abs(trange)*(6.0/20.0); //
  float nhit_thresh   = TMath::Abs(trange)*nHitThreshFact;
  //looop over hits
  for (int i=0;i<nUse;i++){
    if (trange>=0){
      if ((ttc[0][i]>thit)&&(ttc[0][i]<thit+trange)) nhits++;
    }
    else {
      if ((ttc[0][i]<thit)&&(ttc[0][i]>thit+trange)) nhits++;
    }
  }
  if (nhits>=nhit_thresh){
    return 1;
  }
  else {
    return 0;
    cout<<"is NOT a valid hit!!"<<endl;
  }

}

int spliTChan::mergeClusters(){
  if (Cluster_nCand==1) return 0;
  else {
    int mergingClusters = 1;
    int clusterMerged;
    int nClustersMerged = 0;
    while (mergingClusters==1){
      clusterMerged = 0;
      if (Cluster_nCand<=1) break;
      cout<<"looping over "<<Cluster_nCand<<" clusters"<<endl;
      for (int i=0;i<Cluster_nCand-1;i++){
        cout<<"checking "<<i<<"th cluster"<<endl;
        clusterMerged=checkCluster(i);
        cout<<"result: "<<clusterMerged<<endl;
        if (clusterMerged==1){
          nClustersMerged++;
          break;
        } 
      }
      mergingClusters=clusterMerged;
      cout<<"mergingclusters is: "<<mergingClusters<<endl;
      cout<<"c merged is: "<<clusterMerged<<endl;
    } 
    cout<<"merged "<<nClustersMerged<<" clusters!"<<endl;
    return 0;
  }        
}

int spliTChan::checkCluster(int ncluster){
  //check if cluster 2 does NOT overlap with cluster one:
  if (Cluster_tStart[ncluster+1]>Cluster_tEnd[ncluster]){
    return 0;
  }
  else {
    cout<<"merging cluster: "<<ncluster<<endl;
    cout<<"cluster start: "<<Cluster_tStart[ncluster]<<endl;
    cout<<"cluster end  : "<<Cluster_tEnd[ncluster]<<endl;
    cout<<"cluster 2 start: "<<Cluster_tStart[ncluster+1]<<endl;
    cout<<"cluster 2 end: "<<Cluster_tEnd[ncluster+1]<<endl;

    //merge clusters
    Cluster_tStart[ncluster] = TMath::Min(Cluster_tStart[ncluster],Cluster_tStart[ncluster+1]);
    Cluster_tEnd[ncluster] = TMath::Max(Cluster_tEnd[ncluster],Cluster_tEnd[ncluster+1]);
    this->RemoveCluster(ncluster+1);
//    Cluster_nCand--;
    return 1;
  }
}

//amis
//Peak-to-cluster (P2C) cluster finding algorithm
void spliTChan::MakeClusters(int npeaks, float* peak_t0, double* vertex, int nevt){
  //set up fixed time window//
 // cout<<"called MakeClusters"<<endl;
 // float window_start = 180;  //nominal value 180
 // float window_end   = 800;   //nominal value 800
  float window_start = fixedTBeforePk;
  float window_end   = fixedTAfterPk;
  ////
  //?//
  float trange1   = 0.0;
  float trange2   = 0.0;
  //
  //half width of window centered at t0 to count hits
 // cout<<"0th peak time is: "<<peak_t0[0];
  Cluster_nCand = npeaks; //set # of clusters equal to # peaks
  //
  float htmax  = window_end + peak_t0[npeaks-1]; //set maximum possible tof hit time
  float htmin  = peak_t0[0]-window_start;  //set minimum possible tof hit time
  float tmin[npeaks]; //array of window start times
  float tmax[npeaks]; //array of window end times
 // cout<<"tmin, tmax: "<<tmin<<","<<tmax<<endl;
 // cout<<"NEVENT IS: "<<nevt<<endl;
  GetHitInfo(vertex);
//  cout<<"got hit info"<<endl;
  //make sure ttc array has been filled:
  if (ttc[0][0] == 0) return;
  if (ttc[1][0] == 0) return;
 // cout<<"nozero first elements in ttc"<<endl;
 // cout<<"ttc[0][0]: "<<ttc[0][0]<<endl;
//  cout<<"ttc[1][1]: "<<ttc[1][1]<<endl;
  double thit;
  float c_start_time;
  float c_end_time;
  float tof_start_time;
  float tof_end_time;

  for (int j=0;j<Cluster_nCand;j++){
    //loop over hits in peak to find cluster start, end times
    //define window around peak from which to accept hits:
    tmin[j] = peak_t0[j] - window_start;
    tmax[j] = peak_t0[j] + window_end;
    c_start_time = 1000000.;
    c_end_time   = 0.;
  //  cout<<"making cluster from peak "<<Cluster_nCand<<endl;
    for (int i=0;i<nUse;i++){
 //       cout<<"checking hit "<<i<<endl;
 //     cout<<"raw: "<<ttc[0][i]<<endl;
 //     cout<<"cor: "<<ttc[1][i]<<endl;
      if (ttc[1][i]<tmin[j]) continue;                 //skip hits outside time window?
      if (ttc[1][i]>tmax[j]) continue;
      if (ttc[0][i]<c_start_time){
 //       cout<<"calling isValidHit"<<endl;
        if (isValidHit(ttc[0][i],trange1)==1){
          c_start_time = ttc[0][i];
          tof_start_time = ttc[1][i];
        }
      }
      if (ttc[0][i]>c_end_time){
        if (isValidHit(ttc[0][i],trange2)==1){
          c_end_time = ttc[0][i];
          tof_end_time = ttc[1][i];
        }
      }
    }
    // make sure time window is valid
    if (c_start_time>=c_end_time){
      //repeat algorithm without valid checking
      for (int k=0;k<nUse;k++){
        if (ttc[1][k]<tmin[j]) continue; //skip hits outside time window?
        if (ttc[1][k]>tmax[j]) continue;
        if (ttc[0][k]<c_start_time){
          c_start_time = ttc[k][0];
        }
        if (ttc[0][k]>c_end_time){
          c_end_time = ttc[0][k];
        }
      }
   }
    //set cluster properties
    Cluster_tStart[j] = c_start_time;
    Cluster_tEnd[j]   = c_end_time;
//    tof_cluster_tstart[j] = tof_start_time;
//    tof_cluster_tend[j]   = tof_end_time;
    Cluster_goodflag[j] = 1;
  }
//  cout<<"first cluster: "<<Cluster_tStart[0]<<", "
//      <<Cluster_tEnd[0]<<endl;
  //merge clusters
  mergeClusters();
  //calc cluster properties
  CalcClusterProperties();
  //fill common blocs 
  FillClusterCommon();
}


void spliTChan::SetATMClusters(){
  
  cout << "Setting clusters based on ATM gates..." << endl;
  
  Cluster_nCand = fQEventManager::Get()->n_Gate; // # of ATM gates
  
  for (int j=0;j<Cluster_nCand;j++){
    
    //set cluster properties
    Cluster_tStart[j] = fQEventManager::Get()->t_Gate[j][0];
    Cluster_tEnd[j]   = fQEventManager::Get()->t_Gate[j][1];
    
    Cluster_goodflag[j] = 1;
  }
  
  //calc cluster properties
  CalcClusterProperties();
  //fill common blocs 
  FillClusterCommon();
}