diff --git a/DataModel/RecoCluster.cpp b/DataModel/RecoCluster.cpp
index c63855600..db6693efd 100644
--- a/DataModel/RecoCluster.cpp
+++ b/DataModel/RecoCluster.cpp
@@ -1,4 +1,4 @@
-#include "RecoCluster.h"
+#include "RecoCluster.h"
 #include "ANNIERecoObjectTable.h"
 #include <algorithm>
 
@@ -16,37 +16,425 @@ RecoCluster::~RecoCluster()
 
 void RecoCluster::Reset()
 {
-  int j=fDigitList.size();
-  for(int i=0;i<j;i++){
+  int j=fDigitList->size();
+  /*for (int i = 0; i<j; i++) {
     delete fDigitList[i];
-  }
-  fDigitList.clear();
+  }*/
+  fDigitList->clear();
 }
 
-static bool CompareTimes(RecoDigit *rd1, RecoDigit *rd2)
+static bool CompareTimes(RecoDigit rd1, RecoDigit rd2)
 {
-  return ( rd1->GetCalTime() > rd2->GetCalTime() );
+  return ( rd1.GetCalTime() > rd2.GetCalTime() );
 }
 
 void RecoCluster::SortCluster()
 {
-  sort(fDigitList.begin(), fDigitList.end(), CompareTimes);
+  sort(fDigitList->begin(), fDigitList->end(), CompareTimes);
 }
 
-void RecoCluster::AddDigit(RecoDigit* digit)
+void RecoCluster::AddDigit(RecoDigit digit)
 {
-  fDigitList.push_back(digit);
+  //digit.AddCluster(fClusterMode);
+    fDigitList->push_back(digit);
 }
 
-RecoDigit* RecoCluster::GetDigit(int n)
+RecoDigit RecoCluster::GetDigit(int n)
 {
-  return (RecoDigit*)(fDigitList.at(n));
+  return (RecoDigit)(fDigitList->at(n));
 }
 
 int RecoCluster::GetNDigits()
 {
-  return fDigitList.size();
+  return fDigitList->size();
 }
 
+void RecoCluster::SetClusterMode(int cmode)
+{
+  fClusterMode = cmode;
+}
+
+int RecoCluster::GetClusterMode()
+{
+  return fClusterMode;
+}
+
+bool RecoCluster::CheckFilter() {
+    map<int,double> ModeDigitMap;
+    vector<int> DigitModes;
+    double baseline;
+
+    fIsFiltered=true;
+    for (int i=0;i<fDigitList->size();i++) {
+        RecoDigit adigit=fDigitList->at(i);
+        DigitModes=adigit.GetClusteredModes();
+        for (int j = 0; j < DigitModes.size(); j++) {
+            if(!ModeDigitMap.count(DigitModes.at(j))) ModeDigitMap.emplace(DigitModes.at(j),adigit.GetCalCharge());
+            else ModeDigitMap.at(DigitModes.at(j))+=adigit.GetCalCharge();
+        }
+        
+    }
+    if(ModeDigitMap.size()==0) return fIsFiltered;
+    baseline=ModeDigitMap.at(fClusterMode);
+    for (pair<int, double> apair : ModeDigitMap) {
+        if (apair.first>=fClusterMode) continue;
+        if(apair.second/baseline>0.6) fIsFiltered=false;
+    }
+    return fIsFiltered;
+}
+
+void RecoCluster::CalcTime() {
+	double sum = 0;
+	for (int i=0; i<fDigitList->size(); i++) {
+		sum += fDigitList->at(i).GetCalTime();
+	}
+	clusterTime = sum / GetNDigits();
+
+}
+
+void RecoCluster::CalcCharge() {
+    double sum=0;
+    for (int i = 0; i < fDigitList->size(); i++) {
+        sum += fDigitList->at(i).GetCalCharge();
+    }
+    clusterCharge = sum;
+}
+
+void RecoCluster::CalcCB() {
+    //calculate unfiltered CB
+    double total_Q = 0;
+    double total_QSquared = 0;
+    for (int i = 0; i < fDigitList->size(); i++) {
+        //if(unfilteredDigits->at(i)->GetDigitType()==RecoDigit::PMT8inch){
+        if (fDigitList->at(i).GetFilterStatus()) {
+            double tube_charge = fDigitList->at(i).GetCalCharge();
+            total_Q += tube_charge;
+            total_QSquared += (tube_charge * tube_charge);
+            //}
+        }
+    }
+    //FIXME: Need a method to have the 123 be equal to the number of operating detectors
+    double charge_balance = sqrt((total_QSquared) / (total_Q * total_Q) - (1. / 123.));
+
+    clusterCB=charge_balance;
+}
+
+//Calculate Angular span
+// - AS0 is the maximum angle between any two hit PMTs
+// - AS1 is the largest angle between the greatest hit by charge and any of the other hits
+// - ASC is a charge-weighted version of AS0
+void RecoCluster::CalcAS() {
+
+        double max_angle = 0;
+        double angle;
+        Position i_position, j_position;
+        double max_charge = 0;
+        double max_angle2 = 0;
+        int max_index = 0;
+        Position max_position;
+        double max_chargeangle = 0;
+        double chargeangle;
+        double iq, jq;
+
+        for (int i = 0; i < fDigitList->size(); i++) {
+            i_position = fDigitList->at(i).GetPosition();
+            for (int j = 0; j < fDigitList->size(); j++) {
+                if (j == i)continue;
+                j_position = fDigitList->at(j).GetPosition();
+
+                angle = j_position.Angle(i_position);
+                if (angle > max_angle)max_angle = angle;
+
+            }
+            
+            if (fDigitList->at(i).GetCalCharge() > max_charge) {
+                max_charge = fDigitList->at(i).GetCalCharge();
+                max_index = i;
+                max_position = fDigitList->at(i).GetPosition();
+            }
+
+            iq = fDigitList->at(i).GetCalCharge();
+            for (int j = 0; j < fDigitList->size(); j++) {
+                if (j == i)continue;
+                j_position = fDigitList->at(j).GetPosition();
+                jq = fDigitList->at(j).GetCalCharge();
+                chargeangle = iq * jq * j_position.Angle(i_position);
+                if (chargeangle > max_chargeangle)max_chargeangle = chargeangle;
+            }
+        }
+        AS0 = max_angle;
+
+        for (int i = 0; i < fDigitList->size(); i++) {
+            if (i == max_index)continue;
+            i_position = fDigitList->at(i).GetPosition();
+            angle = i_position.Angle(max_position);
+            if (angle > max_angle2)max_angle2 = angle;
+        }
+        AS1 = max_angle2;
+
+        ASC=max_chargeangle;
+
+}
+
+//Retrieve Angular Span value
+double RecoCluster::GetAS(int mode) {
+    if(mode==0)return AS0;
+    else if(mode==1)return AS1;
+    else return 0;
+}
+
+void RecoCluster::CalcParameters() {
+    this->CalcTime();
+    this->CalcCharge();
+    this->CalcCB();
+    this->CalcAS();
+    this->CalcSA();
+    this->CalcAW();
+    this->CalcCV();
+    this->CalcAMD();
+    this->CalcPlanaritySphericity();
+    this->CalcTR();
+
+    //Other needed parameters here
+}
+
+void RecoCluster::SetParticle(int pID, int pPDG, double eff, double pur) {
+    bestParticleID=pID;
+    bestParticlePDG=pPDG;
+    efficiency=eff;
+    purity=pur;
+}
+
+int RecoCluster::calcBestParent() {
+    std::map<int,double> ParticletoClusterCharge;
+    int parentCandidate;
+    int maxIndex=0;
+    double digitCharge;
+    double maxCharge=0;
+    int tempPDG = -5;
+    for (int i = 0; i < fDigitList->size(); i++) {
+        RecoDigit i_digit=fDigitList->at(i);
+        for(int j=0; j<i_digit.GetParents().size(); j++){
+        parentCandidate=i_digit.GetParents().at(i);
+        digitCharge=i_digit.GetCalCharge();
+        if (ParticletoClusterCharge.count(parentCandidate) == 0) {
+            ParticletoClusterCharge.emplace(parentCandidate,digitCharge);
+        }
+        else {
+            ParticletoClusterCharge.at(parentCandidate)+=digitCharge;
+        }
+        }
+    }
+    for (std::pair<int, double>apair : ParticletoClusterCharge) {
+        if (apair.second > maxCharge) {
+            maxCharge=apair.second;
+            maxIndex=apair.first;
+            
+        }
+    }
 
+    //cout << "Particle " << maxIndex << " wins with charge " << maxCharge << endl;
+
+    bestParticleID=maxIndex;
+    purity=maxCharge/clusterCharge;
+
+    return maxIndex;
+}
+
+void RecoCluster::CalcCV() {
+    Position ud;
+    ChargeVector.SetX(0);
+    ChargeVector.SetY(0);
+    ChargeVector.SetZ(0);
+    for (int i = 0; i < fDigitList->size(); i++) {
+        for (int j = i + 1; j < fDigitList->size(); j++) {
+            ud=(fDigitList->at(i).GetPosition() - fDigitList->at(j).GetPosition()).Unit();
+            ChargeVector+=(fDigitList->at(i).GetCalCharge()*fDigitList->at(j).GetCalCharge())*ud;
+        }
+    }
+}
+
+void RecoCluster::CalcAMD() {
+    double min_distance = 9999;
+    double ave_min_distance=0;
+    double distance;
+    for (int i = 0; i < fDigitList->size(); i++) {
+        Position ipos=fDigitList->at(i).GetPosition();
+        for (int j = 0; j < fDigitList->size(); j++) {
+            if(j == i) continue;
+            Position jpos=fDigitList->at(j).GetPosition();
+            distance = (jpos-ipos).Mag();
+            if(distance<min_distance) min_distance=distance;
+
+        }
+        ave_min_distance+=min_distance;
+        min_distance=9999;
+    }
+    ave_min_distance/=fDigitList->size();
+    AMD=ave_min_distance;
+}
+
+//spatial average: Charge-averaged location of all 
+void RecoCluster::CalcSA() {
+   double aveX=0;
+   double aveY=0;
+   double aveZ=0;
+   for (int i = 0; i < fDigitList->size(); i++) {
+       Position ipos=fDigitList->at(i).GetPosition();
+       double iq=fDigitList->at(i).GetCalCharge();
+       aveX += ipos.X() * iq;
+       aveY += ipos.Y() * iq;
+       aveZ += ipos.Z() * iq;
+   }
+   aveX/=clusterCharge;
+   aveY/=clusterCharge;
+   aveZ/=clusterCharge;
+   SpatialAverage=Position(aveX,aveY,aveZ);
+}
+
+//Containing Angle: angle from SpatialAverage unit vector which contains 68% of the total charge of the cluster
+void RecoCluster::CalcAW() {        //AW: Angular Width?  1-sigma angular width?
+    std::map<double,double> angle_charge_map;
+    double angle,charge;
+    double contained=0;
+    for (int i = 0; i < fDigitList->size(); i++) {
+        RecoDigit idigit=fDigitList->at(i);
+        angle=SpatialAverage.Angle(idigit.GetPosition());
+        charge = idigit.GetCalCharge();
+        if(angle_charge_map.count(angle)==0)
+            angle_charge_map.emplace(angle,charge);
+        else
+            angle_charge_map.at(angle)+=charge;
+    }
+    for (std::pair<double,double> apair: angle_charge_map) {
+        AngularWidth=apair.first;
+        contained+=apair.second/clusterCharge;
+        if(contained>0.68)break;
+    }
+    return;
+}
+
+//Planarity and Sphericity: measure of how planar or spherical cluster is; currently in-process of translating from reference
+void RecoCluster::CalcPlanaritySphericity()
+{
+    double P = 0.0; double S = 0.0;
+
+    // 1) Charge-weighted centroid
+    Position mean(0, 0, 0);
+    double wsum = 0.0;
+    for (int i = 0; i<fDigitList->size(); i++) {
+        double w = fDigitList->at(i).GetCalCharge();
+        mean += w * fDigitList->at(i).GetPosition();
+        wsum += w;
+    }
+    if (wsum <= 0.0) return;
+    mean *= (1.0 / wsum);
+
+    // 2) Charge-weighted covariance (symmetric 3x3)
+    //    M = (1/wsum) * sum_i w_i * (r_i)(r_i)^T, with r_i = pos_i - mean
+    TMatrixDSym cov(3);
+    cov.Zero();
+    for (int i = 0; i < fDigitList->size(); i++) {
+        RecoDigit d = fDigitList->at(i);
+        double w = d.GetCalCharge();
+        Position r = d.GetPosition() - mean;
+        // Outer product r*r^T scaled by w
+        cov(0, 0) += w * r.X() * r.X();
+        cov(0, 1) += w * r.X() * r.Y();
+        cov(0, 2) += w * r.X() * r.Z();
+        cov(1, 1) += w * r.Y() * r.Y();
+        cov(1, 2) += w * r.Y() * r.Z();
+        cov(2, 2) += w * r.Z() * r.Z();
+    }
+    cov(1, 0) = cov(0, 1);
+    cov(2, 0) = cov(0, 2);
+    cov(2, 1) = cov(1, 2);
+    cov *= (1.0 / wsum);
+
+    // 3) Eigen-decomposition (symmetric → real eigen system)
+    TMatrixDSymEigen eig(cov);
+    TVectorD evals = eig.GetEigenValues(); // typically ascending, but we’ll sort to be safe
+
+    // Copy to std::vector and sort descending
+    std::vector<double> lambda = { evals[0], evals[1], evals[2] };
+    std::sort(lambda.begin(), lambda.end(), std::greater<double>());
+
+    const double lambda1 = lambda[0]; // largest
+    const double lambda2 = lambda[1];
+    const double lambda3 = lambda[2]; // smallest
+    const double sum = lambda1 + lambda2 + lambda3;
+
+    if (sum <= 0.0) return;
+
+    // 4) Planarity and Sphericity
+    //    P = (lambda2 - lambda3) / (lambda1 + lambda2 + lambda3)
+    //    S = (3/2) * lambda3 / (lambda1 + lambda2 + lambda3)
+    P = (lambda2 - lambda3) / sum;
+    S = 1.5 * (lambda3 / sum);
+
+    Planarity=P;
+    Sphericity=S;
+}
+
+double RecoCluster::CalcBeta(int order) {
+    double Beta=0;
+    Position ipos,jpos;
+    double angle;
+    int NDigits=fDigitList->size();
+    TF1 LegendreFunction("legendrefunction","ROOT::Math::legendre([0],x)",-1,1);
+    LegendreFunction.SetParameter(0,order);
+    for (int i = 0; i < fDigitList->size(); i++) {
+        for (int j = i + 1; j < fDigitList->size(); j++) {
+            ipos=fDigitList->at(i).GetPosition();
+            jpos=fDigitList->at(j).GetPosition();
+            angle=ipos.Angle(jpos);
+            Beta+=2 * LegendreFunction.Eval(cos(angle)) / (NDigits*(NDigits-1));
+        }
+    }
+
+    if(BetaParameters.count(order)) BetaParameters.at(order)=Beta;
+    else BetaParameters.emplace(order, Beta);
+
+    return Beta;
+}
+
+double RecoCluster::GetBeta(int order) {
+    if(BetaParameters.count(order))return BetaParameters.at(order);
+    else return -999;
+}
+
+void RecoCluster::CalcTR() {
+    double firsttime=fDigitList->at(0).GetCalTime();
+    double lasttime=fDigitList->at(fDigitList->size()-1).GetCalTime();
+    double maxtime=fDigitList->at(0).GetCalTime();
+    double maxcharge=fDigitList->at(0).GetCalCharge();
+    
+    for (int i = 0; i < fDigitList->size(); i++) {
+        RecoDigit adigit=fDigitList->at(i);
+        if(adigit.GetCalTime()<firsttime)firsttime=adigit.GetCalTime();
+        if(adigit.GetCalTime()>lasttime)lasttime=adigit.GetCalTime();
+        if (adigit.GetCalCharge() > maxcharge) {
+            maxtime=adigit.GetCalTime();
+            maxcharge=adigit.GetCalCharge();
+        }
+    }
+
+    TRTotal=lasttime-firsttime;
+    if(clusterTime-firsttime>lasttime-clusterTime)TRCluster=clusterTime-firsttime;
+    else TRCluster=lasttime-clusterTime;
+    if(maxtime-firsttime>lasttime-maxtime)TRQ=maxtime-firsttime;
+    else TRQ=lasttime-maxtime;
+
+    TRRatioTQ=TRTotal/TRQ;
+    TRRatioTC=TRTotal/TRCluster;
+    TRRatioQC=TRQ/TRCluster;
+}
 
+/*void RecoCluster::CleanDigits() {
+    for (int i = 0; i < fDigitList->size(); i++) {
+        if (fDigitList->at(i) != nullptr) {
+            delete fDigitList->at(i);
+            fDigitList->at(i)=nullptr;
+        }
+    }
+}*/
\ No newline at end of file
diff --git a/DataModel/RecoCluster.h b/DataModel/RecoCluster.h
index 4afb14268..d338c392a 100644
--- a/DataModel/RecoCluster.h
+++ b/DataModel/RecoCluster.h
@@ -3,7 +3,15 @@
 #ifndef RECOCLUSTER_H
 #define RECOCLUSTER_H
 #include "RecoDigit.h"
+#include "Position.h"
 #include<SerialisableObject.h>
+#include "TVector3.h"
+#include "TMatrixDSym.h"
+#include "TMatrixDSymEigen.h"
+#include "TVectorD.h"
+#include "TMath.h"
+#include "TF1.h"
+#include "Math/IFunction.h"
 
 class RecoCluster : public SerialisableObject {
 
@@ -16,19 +24,96 @@ class RecoCluster : public SerialisableObject {
   void Reset();
   void SortCluster();
 
-  void AddDigit(RecoDigit* digit);
+  void AddDigit(RecoDigit digit);
+  inline void SetDigits(vector<RecoDigit>* indigits){fDigitList=indigits;}
 
-  RecoDigit* GetDigit(int n);
+  RecoDigit GetDigit(int n);
+  
+  void SetClusterMode(int cmode);
+  
+  int GetClusterMode();
+  
+  std::vector<RecoDigit>* GetDigitList() {return this->fDigitList;}
+  
   int GetNDigits();
   
   bool Print() {
 		cout<<"Number of digits in this cluster : "<<GetNDigits()<<endl;
+		cout<<"Clustering mode : "<<GetNDigits()<<endl;
 		return true;
 	}
 
+  inline double GetTime() { return clusterTime; }
+  inline double GetCharge() { return clusterCharge; }
+
+  void CalcTime();
+  void CalcCharge();
+
+  inline double GetCB() { return clusterCB; }
+  inline Position GetCV() {return ChargeVector; }
+
+  void CalcCB();
+  void CalcAS();
+  void CalcCV();
+  void CalcAMD();
+  void CalcSA();
+  void CalcAW();
+  void CalcPlanaritySphericity();
+  double CalcBeta(int order);
+  void CalcTR();
+  void CalcParameters();
+  int calcBestParent();
+  double GetAS(int mode);
+  inline double GetASC(){return ASC;}
+  inline double GetAMD(){return AMD;}
+  inline Position GetSA(){return SpatialAverage;}
+  inline double GetAW(){return AngularWidth;}
+  inline double GetPlanarity(){return Planarity;}
+  inline double GetSphericity(){return Sphericity;}
+  double GetBeta(int order);
+  inline double GetTimeRangeT() {return TRTotal;}
+  inline double GetTimeRangeQ() {return TRQ;}
+  inline double GetTimeRangeC() {return TRCluster;}
+  inline double GetTRRTQ(){return TRRatioTQ;}
+  inline double GetTRRTC() { return TRRatioTC; }
+  inline double GetTRRQC() { return TRRatioQC; }
+
+  void SetParticle(int pID,int pPDG, double eff, double pur);
+  inline int GetBestParent(){return bestParticleID;}
+  inline void SetPDG(int pPDG) {bestParticlePDG=pPDG;}
+  inline int GetPDG(){return bestParticlePDG;}
+  inline double Efficiency(){return efficiency;}
+  inline double Purity(){return purity;}
+
+  bool CheckFilter();
+  inline bool GetFilterStatus(){return fIsFiltered;}
+  //void CleanDigits();
+
   private:
+	  double clusterTime;
+	  double clusterCharge;
+	  double clusterCB;
+	  double AS0, AS1, ASC;
+	  double AMD;
+	  Position ChargeVector;
+	  Position SpatialAverage;
+	  double AngularWidth;
+	  double Planarity=-999;
+	  double Sphericity=-999;
+	  std::map<int,double> BetaParameters;
+	  double TRTotal, TRQ, TRCluster;
+	  double TRRatioTQ, TRRatioTC,TRRatioQC;
+	  int bestParticleID;
+	  int bestParticlePDG;
+	  double efficiency;
+	  double purity;
+	  bool fIsFiltered=0;
+
+  
+  int fClusterMode = -999;
+  std::vector<RecoDigit>* fDigitList;
+  Position TwoDCenter;
 
-  std::vector<RecoDigit*> fDigitList;
   	
   protected:
   template<class Archive> void serialize(Archive & ar, const unsigned int version){
diff --git a/DataModel/RecoDigit.h b/DataModel/RecoDigit.h
index 13d12c4d0..f766d0e45 100644
--- a/DataModel/RecoDigit.h
+++ b/DataModel/RecoDigit.h
@@ -33,15 +33,32 @@ class RecoDigit : public SerialisableObject{
 			fIsFiltered = 1;
 			ANNIERecoObjectTable::Instance()->NewDigit();
 		}
+	RecoDigit(RecoDigit* origin) {
+		serialise=true;
+		fRegion=origin->GetRegion();
+		fPosition = origin->GetPosition();
+		fCalTime = origin->GetCalTime();
+		fCalCharge=origin->GetCalCharge();
+		fDigitType=origin->GetDigitType();
+		fDetectorID=origin->GetDetectorID();
+		fIsFiltered=origin->GetFilterStatus();
+
+	}
 	~RecoDigit() {/*ANNIERecoObjectTable::Instance()->DeleteDigit();*/}
 
+	bool operator<(const RecoDigit other) const {
+		return (fPosition.X() == other.GetPosition().X()) ? fPosition.Y() < other.GetPosition().Y() : fPosition.X() < other.GetPosition().X();
+	}
+
 	inline int                 GetRegion() const {return fRegion;}
 	inline Position						 GetPosition() const {return fPosition;}
 	inline double              GetCalTime() const {return fCalTime;}
 	inline double              GetCalCharge() const {return fCalCharge;}
 	inline bool                GetFilterStatus() const {return fIsFiltered;}
+	inline vector<int>		GetClusteredModes() { return fClusterMode; }
 	inline int         				 GetDigitType() const {return fDigitType;}
 	inline int         				 GetDetectorID() const {return fDetectorID;}
+	inline std::vector<int> GetParents() const {return Parents;}
 	
 	inline void                SetRegion(int reg) {fRegion = reg;}
 	inline void                SetPosition(Position pos){fPosition = pos;}
@@ -52,6 +69,9 @@ class RecoDigit : public SerialisableObject{
 	inline void                SetFilter(bool pass = 1) { fIsFiltered = pass;}
   inline void                ResetFilter() {SetFilter(0);}
   inline void                PassFilter() {SetFilter(1);}
+  inline void				AddCluster(int InClusterMode) {fClusterMode.push_back(InClusterMode); }
+  inline void				UnCluster(){fClusterMode.pop_back(); }
+  inline void SetParents(std::vector<int> parentsin) { Parents = parentsin; }
 
 	bool Print() {
 		cout<<"Region : "<<fRegion<<endl;
@@ -70,8 +90,10 @@ class RecoDigit : public SerialisableObject{
   double fCalTime;
   double fCalCharge; 
   bool fIsFiltered;
+  vector<int> fClusterMode;
   int fDigitType;
   int fDetectorID;
+  std::vector<int> Parents;
 
 	template<class Archive> void serialize(Archive & ar, const unsigned int version){
 		if(serialise){
diff --git a/UserTools/ClusterSearcher/ClusterSearcher.cpp b/UserTools/ClusterSearcher/ClusterSearcher.cpp
new file mode 100644
index 000000000..a882513fa
--- /dev/null
+++ b/UserTools/ClusterSearcher/ClusterSearcher.cpp
@@ -0,0 +1,729 @@
+#include "ClusterSearcher.h"
+
+static ClusterSearcher* fgClusterSearcher = 0;
+
+ClusterSearcher* ClusterSearcher::Instance()
+{
+  if( !fgClusterSearcher ){
+    fgClusterSearcher = new ClusterSearcher();
+  }
+
+  return fgClusterSearcher;
+}
+
+ClusterSearcher::ClusterSearcher():Tool(){}
+	
+ClusterSearcher::~ClusterSearcher() {
+  
+}
+
+bool ClusterSearcher::Initialise(std::string configfile, DataModel &data){
+    if(verbosity)cout<<"Initializing ClusterSearcher"<<endl;
+  /////////////////// Usefull header ///////////////////////
+  if(configfile!="")  m_variables.Initialise(configfile); //loading config file
+  //m_variables.Print();
+
+  m_data= &data; //assigning transient data pointer
+  /////////////////////////////////////////////////////////////////
+  
+
+  
+  /// Get the Tool configuration variables
+  m_variables.Get("verbosity",verbosity);
+  m_variables.Get("IsMC",fisMC);
+  m_variables.Get("ClusterMode",fClusterMode);
+  m_variables.Get("PmtMinPulseHeight", fPmtMinPulseHeight);
+  m_variables.Get("PmtNeighbourRadius", fPmtNeighbourRadius);
+  m_variables.Get("PmtMinNeighbourDigits", fPmtMinNeighbourDigits);
+  m_variables.Get("PmtClusterRadius", fPmtClusterRadius);
+  m_variables.Get("PmtTimeWindowN", fPmtTimeWindowN);	
+  m_variables.Get("PmtTimeWindowC", fPmtTimeWindowC);
+  m_variables.Get("PmtMinHitsPerCluster", fPmtMinHitsPerCluster);	
+  m_variables.Get("LappdMinPulseHeight", fLappdMinPulseHeight);
+  m_variables.Get("LappdNeighbourRadius", fLappdNeighbourRadius);
+  m_variables.Get("LappdMinNeighbourDigits", fLappdMinNeighbourDigits);
+  m_variables.Get("LappdClusterRadius", fLappdClusterRadius );
+  m_variables.Get("LappdTimeWindowN", fLappdTimeWindowN );
+  m_variables.Get("LappdTimeWindowC", fLappdTimeWindowC );
+  m_variables.Get("LappdMinHitsPerCluster", fLappdMinHitsPerCluster );
+  m_variables.Get("MinClusterDigits", fMinClusterDigits );
+  m_variables.Get("SinglePEGains",singlePEgains);
+  m_variables.Get("FirstRun",fFirstRun);
+
+  /// Fill map with settings of ClusterSearcher
+  fClusteringParam = new std::map<std::string,double>;
+  fClusteringParam->emplace("ClusterMode",fClusterMode);
+  fClusteringParam->emplace("PmtMinPulseHeight",fPmtMinPulseHeight);
+  fClusteringParam->emplace("PmtNeighbourRadius",fPmtNeighbourRadius);
+  fClusteringParam->emplace("PmtMinNeighbourDigits",fPmtMinNeighbourDigits);
+  fClusteringParam->emplace("PmtClusterRadius",fPmtClusterRadius);
+  fClusteringParam->emplace("PmtTimeWindowN",fPmtTimeWindowN);
+  fClusteringParam->emplace("PmtTimeWindowC",fPmtTimeWindowC);
+  fClusteringParam->emplace("LappdMinPulseHeight",fLappdMinPulseHeight);
+  fClusteringParam->emplace("LappdNeighbourRadius",fLappdNeighbourRadius);
+  fClusteringParam->emplace("LappdMinNeighbourDigits",fLappdMinNeighbourDigits);
+  fClusteringParam->emplace("LappdClusterRadius",fLappdClusterRadius);
+  fClusteringParam->emplace("LappdTimeWindowN",fLappdTimeWindowN);
+  fClusteringParam->emplace("LappdTimeWindowC",fLappdTimeWindowC);
+  fClusteringParam->emplace("MinClusterDigits",fMinClusterDigits);
+
+  Log("ClusterSearcher " + to_string(fClusterMode) + " configs Initialized", v_debug, verbosity);
+  if (!fisMC){
+    ifstream file_singlepe(singlePEgains.c_str());
+    unsigned long temp_chankey;
+    double temp_gain;
+    while (!file_singlepe.eof()){
+      file_singlepe >> temp_chankey >> temp_gain;
+      if (file_singlepe.eof()) break;
+      pmt_gains.emplace(temp_chankey,temp_gain);
+    }
+    file_singlepe.close();
+    m_data->CStore.Get("pmt_tubeid_to_channelkey",pmt_tubeid_to_channelkey);
+  }
+  Log("ClusterSearcher " + to_string(fClusterMode) + " !MC Initialized", v_debug, verbosity);
+
+  // vector of selected digits
+  fSelectAll = new std::vector<RecoDigit*>;
+  fSelectByPulseHeight = new std::vector<RecoDigit*>;
+  fSelectByNeighbours = new std::vector<RecoDigit*>;
+  fSelectByClusters = new std::vector<RecoDigit*>;
+  // only for test 
+  fSelectByTruthInfo = new std::vector<RecoDigit*>; 
+  // vector of clusters
+  fClusterList = new std::vector<RecoCluster>;
+  fRecoClusters = new std::vector<RecoCluster>;
+
+  //Set clustering parameters in the RecoEvent store
+  std::map<std::string, double>* pre_ClusteringParam = nullptr; // read existing container for parameters
+  bool cluster_parameter_status = m_data->Stores.at("RecoEvent")->Get("ClusteringParameters", pre_ClusteringParam);
+  if(!cluster_parameter_status) m_data->Stores.at("RecoEvent")->Set("ClusteringParameters", fClusteringParam);
+  else {
+    pre_ClusteringParam->insert(fClusteringParam->begin(), fClusteringParam->end());
+    //m_data->Stores.at("RecoEvent")->Set("ClusteringParameters", pre_ClusteringParam);   
+  }
+  Log("ClusterSearcher "+to_string(fClusterMode)+" Initialized",v_debug,verbosity);
+  return true;
+}
+
+
+bool ClusterSearcher::Execute(){
+    if (fRecoClusters) Log("Prelim check1: fRecoClusters size: " + to_string(fRecoClusters->size()), v_debug, verbosity);
+    if (fRecoClusters->size()!=0) fRecoClusters->clear();
+    
+    if(fFirstRun && pre_RecoClusters) pre_RecoClusters->clear();
+  
+  std::string name = "ClusterSearcher::Execute()";
+  Log(name + ": Executing",v_error,verbosity);
+	
+  // print Selecting parameters
+  if(verbosity>v_message) this->PrintParameters();
+	
+  // see if "ANNIEEvent" exists
+  auto get_annieevent = m_data->Stores.count("ANNIEEvent");
+  if(!get_annieevent){
+    Log(name + ": No ANNIEEvent store!",v_error,verbosity); 
+    return false;
+  };
+	
+  /// see if "RecoEvent" exists
+  auto get_recoevent = m_data->Stores.count("RecoEvent");
+  if(!get_recoevent){
+    Log(name + ": No RecoEvent store!",v_error,verbosity); 
+    return false;
+  };
+
+  if (fisMC) {
+    // get true vertex
+    auto get_truevtx = m_data->Stores.at("RecoEvent")->Get("TrueVertex", fTrueVertex);
+    if(!get_truevtx){ 
+      Log(name + ": Error retrieving TrueVertex from RecoEvent!",v_error,verbosity); 
+      return false; 
+    }
+  }
+ 
+  // get digit list
+  auto get_recodigit = m_data->Stores.at("RecoEvent")->Get("RecoDigit",fDigitList);  ///> Get digits from "RecoEvent" 
+  if(!get_recodigit){ 
+    Log("VtxSeedGenerator  Tool: Error retrieving RecoDigits,no digit from the RecoEvent store!",v_error,verbosity); 
+    return false;
+  }
+  // copy to a new digit list
+  std::vector<RecoDigit*>* digits = new std::vector<RecoDigit*>;
+  for(int i=0;i<int(fDigitList->size());i++) {
+  	RecoDigit* recodigitptr = &(fDigitList->at(i));
+    digits->push_back((RecoDigit*)recodigitptr);
+  }
+  
+  // Reset Digit Filter
+  // ============
+  for(int n=0; n<int(digits->size()); n++ ) {
+    digits->at(n)->ResetFilter();
+  }
+  
+  // Set Digit Filter
+  // ==========
+  SelectDigits(digits);
+  
+  fRecoClusters = this->RecoClusters(digits);
+  // Digit clustering done!
+  // =====
+  
+  pre_RecoClusters = nullptr;   // to read existing clusters
+  bool cluster_status = m_data->Stores.at("RecoEvent")->Get("RecoClusters", pre_RecoClusters);
+  
+  Log("ClusterSearcher Tool: Final count of clusters: "+to_string(fRecoClusters->size()),v_debug,verbosity);
+  if(!cluster_status) {
+      pre_RecoClusters=new std::vector<RecoCluster>;
+      for (int i = 0; i < fRecoClusters->size(); i++) {
+          pre_RecoClusters->push_back(fRecoClusters->at(i));
+      }
+    fIsHitClusteringDone = true;
+    Log("Hit Clustering done",v_debug,verbosity);
+    m_data->Stores.at("RecoEvent")->Set("HitClusteringDone", fIsHitClusteringDone);
+    m_data->Stores.at("RecoEvent")->Set("RecoClusters", pre_RecoClusters); 
+  }
+  else {
+    pre_RecoClusters->insert(pre_RecoClusters->end(), fRecoClusters->cbegin(), fRecoClusters->cend()); // append new clusters to the existing cluster list
+    //m_data->Stores.at("RecoEvent")->Set("RecoClusters", pre_RecoClusters); // save updated clusters
+    //fRecoClusters->clear();
+  }
+
+  if (!fFirstRun && pre_RecoClusters->size() == 0) {
+      Log("ClusterSearcher: NO CLUSTERS FOUND IN EVENT!",v_debug,verbosity);
+      Log("Digit list size: "+to_string(fDigitList->size()),v_debug,verbosity);
+  }
+
+  delete digits; digits = 0;
+  return true;
+}
+
+bool ClusterSearcher::Finalise(){
+  //delete fClusteringParam; fClusteringParam = 0;      //Will be deleted by the store, don't manually delete
+  delete fSelectAll; fSelectAll = 0;
+  delete fSelectByPulseHeight; fSelectByPulseHeight = 0;
+  delete fSelectByNeighbours; fSelectByNeighbours = 0;
+  delete fSelectByClusters; fSelectByClusters = 0;
+  //delete fRecoClusters; fRecoClusters = 0;    //Will be deleted by the store, don't manually delete
+  delete fClusterList; fClusterList = 0;
+  // for test
+  delete fSelectByTruthInfo; fSelectByTruthInfo = 0;
+  return true;
+}
+
+/*void ClusterSearcher::Config(int config)
+{
+  ClusterSearcher::Instance()->SetConfig(config);
+}*/
+
+void ClusterSearcher::PmtMinPulseHeight(double min)
+{
+  ClusterSearcher::Instance()->SetPmtMinPulseHeight(min);
+}
+  
+void ClusterSearcher::PmtNeighbourRadius(double radius)
+{
+  ClusterSearcher::Instance()->SetPmtNeighbourRadius(radius);
+}
+  
+void ClusterSearcher::PmtNeighbourDigits(int digits)
+{
+  ClusterSearcher::Instance()->SetPmtNeighbourDigits(digits);
+}
+  
+void ClusterSearcher::PmtClusterRadius(double radius)
+{
+  ClusterSearcher::Instance()->SetPmtClusterRadius(radius);
+}
+  
+void ClusterSearcher::MinClusterDigits(int digits)
+{
+  ClusterSearcher::Instance()->SetMinClusterDigits(digits);
+}
+  
+void ClusterSearcher::PmtTimeWindowN(double windowN)
+{
+  ClusterSearcher::Instance()->SetPmtTimeWindowNeighbours(windowN);
+}
+
+void ClusterSearcher::PmtTimeWindowC(double windowC)
+{
+  ClusterSearcher::Instance()->SetPmtTimeWindowClusters(windowC);
+}
+
+void ClusterSearcher::PrintParameters()
+{
+  std::cout << " *** ClusterSearcher::PrintParameters() *** " << std::endl;
+  
+  std::cout << "  Clustering Parameters: " << std::endl
+            << "   ClusterMode = " << fClusterMode<< std::endl
+            << "   PmtMinPulseHeight = " << fPmtMinPulseHeight << std::endl
+            << "   PmtNeighbourRadius = " << fPmtNeighbourRadius << std::endl
+            << "   PmtMinNeighbourDigits = " << fPmtMinNeighbourDigits << std::endl
+            << "   PmtClusterRadius = " << fPmtClusterRadius << std::endl
+            << "   PmtTimeWindowN = " << fPmtTimeWindowN << std::endl
+	          << "   PmtTimeWindowC = " << fPmtTimeWindowC << std::endl
+            << "   LappdMinPulseHeight = " << fLappdMinPulseHeight << std::endl
+            << "   LappdNeighbourRadius = " << fLappdNeighbourRadius << std::endl
+            << "   LappdMinNeighbourDigits = " << fLappdMinNeighbourDigits << std::endl
+            << "   LappdClusterRadius = " << fLappdClusterRadius << std::endl
+            << "   LappdTimeWindowN = " << fLappdTimeWindowN << std::endl
+	          << "   LappdTimeWindowC = " << fLappdTimeWindowC << std::endl
+	          << "   MinClusterDigits = " << fMinClusterDigits << std::endl;
+}
+
+void ClusterSearcher::Reset()
+{
+
+  return;
+}
+
+// Set all filter status to 0 (default is 1)
+std::vector<RecoDigit*>* ClusterSearcher::ResetDigits(std::vector<RecoDigit*>* DigitList)
+{
+  for(int idigit=0; idigit<int(DigitList->size()); idigit++ ){
+    RecoDigit* recoDigit = (RecoDigit*)(DigitList->at(idigit));
+    recoDigit->ResetFilter();
+  }
+
+  return DigitList;
+}
+
+std::vector<RecoDigit*>* ClusterSearcher::SelectDigits(std::vector<RecoDigit*>* DigitList)
+{
+  for(int idigit=0; idigit<int(DigitList->size()); idigit++ ){
+    RecoDigit* recoDigit = (RecoDigit*)(DigitList->at(idigit));
+    recoDigit->PassFilter();
+  }
+
+  return DigitList;
+}
+
+void ClusterSearcher::SelectDigits(std::vector<RecoDigit>* DigitList)
+{
+    for (int idigit = 0; idigit<int(DigitList->size()); idigit++) {
+        RecoDigit recoDigit = (RecoDigit)(DigitList->at(idigit));
+        recoDigit.PassFilter();
+    }
+
+    return;
+}
+
+std::vector<RecoDigit*>* ClusterSearcher::SelectAll(std::vector<RecoDigit*>* DigitList)
+{
+	std::string name = "ClusterSearcher::SelectAll() ";
+  // clear vector of selected digits
+  // ==============================
+  fSelectAll->clear();
+  // select all digits
+  // =================
+  for(int idigit=0; idigit<int(DigitList->size()); idigit++ ){
+    RecoDigit* recoDigit = (RecoDigit*)(DigitList->at(idigit));
+    fSelectAll->push_back(recoDigit);
+  }
+
+  // return vector of selected digits
+  // ================================
+  if(verbosity>v_message) std::cout << name << "  select all: " << fSelectAll->size() << std::endl;
+  
+  return fSelectAll;
+}
+
+std::vector<RecoDigit*>* ClusterSearcher::SelectByPulseHeight(std::vector<RecoDigit*>* DigitList)
+{
+  std::string name = "ClusterSearcher::SelectByPulseHeight() ";
+  // clear vector of selected digits
+  // ===============================
+  fSelectByPulseHeight->clear();
+
+  // Select by pulse height
+  // ======================
+  for(int idigit=0; idigit<int(DigitList->size()); idigit++ ){
+    RecoDigit* recoDigit = (RecoDigit*)(DigitList->at(idigit));
+    double qep = recoDigit->GetCalCharge();
+    if (!fisMC){
+      int pmtid = recoDigit->GetDetectorID();
+      unsigned long chankey = pmt_tubeid_to_channelkey[pmtid];
+      if (pmt_gains[chankey]>0) qep/=pmt_gains[chankey];
+    }
+    int detType = recoDigit->GetDigitType();
+    if(detType == RecoDigit::lappd_v0) {
+    	if( qep>fLappdMinPulseHeight ){
+        fSelectByPulseHeight->push_back(recoDigit);
+      }
+    }
+    else if(detType == RecoDigit::PMT8inch) {
+    	if( qep>fPmtMinPulseHeight ){
+        fSelectByPulseHeight->push_back(recoDigit);
+      }
+    }
+    else Log(name + "<error:> detector type doesn't exist!",v_error,verbosity);
+  }
+
+
+  // return vector of selected digits
+  // ================================
+  if(verbosity>v_message) std::cout <<name << "  Select by pulse height: " << fSelectByPulseHeight->size() << std::endl;
+  
+  return fSelectByPulseHeight;
+}
+
+std::vector<RecoDigit*>* ClusterSearcher::SelectByNeighbours(std::vector<RecoDigit*>* DigitList)
+{
+	std::string name = "ClusterSearcher::SelectByNeighbours() ";
+  // clear vector of Selected digits
+  // ===============================
+  fSelectByNeighbours->clear();
+
+  // create array of neighbours
+  // ==========================
+  int Ndigits = DigitList->size();
+
+  if( Ndigits<=0 ){
+    return fSelectByNeighbours;
+  }
+
+  int* numNeighbours = new int[Ndigits];
+
+  for( int idigit=0; idigit<Ndigits; idigit++ ){
+    numNeighbours[idigit] = 0;
+  }
+
+  // count number of neighbours
+  // ==========================
+  for(int idigit1=0; idigit1<int(DigitList->size()); idigit1++ ){
+  	RecoDigit* fdigit1 = (RecoDigit*)(DigitList->at(idigit1));
+  	TString digit1Type = fdigit1->GetDigitType();
+    for(int idigit2=idigit1+1; idigit2<int(DigitList->size()); idigit2++ ){
+      RecoDigit* fdigit2 = (RecoDigit*)(DigitList->at(idigit2));
+      TString digit2Type = fdigit2->GetDigitType();
+
+      double dx = fdigit1->GetPosition().X() - fdigit2->GetPosition().X();
+      double dy = fdigit1->GetPosition().Y() - fdigit2->GetPosition().Y();
+      double dz = fdigit1->GetPosition().Z() - fdigit2->GetPosition().Z();
+      double dt = fdigit1->GetCalTime() - fdigit2->GetCalTime();
+      double drsq = dx*dx + dy*dy + dz*dz;
+      
+      if(digit1Type == RecoDigit::PMT8inch && digit2Type == RecoDigit::PMT8inch) {
+        if( drsq>0.0
+         && drsq<fPmtNeighbourRadius*fPmtNeighbourRadius
+         && fabs(dt)<fPmtTimeWindowN ){
+          numNeighbours[idigit1]++;
+          numNeighbours[idigit2]++;
+        }	
+      }
+      
+      if(digit1Type == RecoDigit::lappd_v0 && digit2Type == RecoDigit::lappd_v0) {
+        if( drsq>0.0
+         && drsq<fLappdNeighbourRadius*fLappdNeighbourRadius
+         && fabs(dt)<fLappdTimeWindowN ){
+          numNeighbours[idigit1]++;
+          numNeighbours[idigit2]++;
+        }	
+      }
+      
+      if(digit1Type == RecoDigit::PMT8inch && digit2Type == RecoDigit::lappd_v0) {
+        if( drsq>0.0
+            && drsq<fPmtNeighbourRadius*fPmtNeighbourRadius
+            && fabs(dt)<fPmtTimeWindowN ){
+          numNeighbours[idigit1]++;
+        }	
+        if( drsq>0.0
+            && drsq<fLappdNeighbourRadius*fLappdNeighbourRadius
+            && fabs(dt)<fLappdTimeWindowN ){
+          numNeighbours[idigit2]++;
+        }	   
+      }
+
+      if(digit1Type == RecoDigit::lappd_v0 && digit2Type == RecoDigit::PMT8inch) {
+        if( drsq>0.0
+            && drsq<fLappdNeighbourRadius*fLappdNeighbourRadius
+            && fabs(dt)<fLappdTimeWindowN ){
+          numNeighbours[idigit1]++;
+        }	
+        if( drsq>0.0
+            && drsq<fPmtNeighbourRadius*fPmtNeighbourRadius
+            && fabs(dt)<fPmtTimeWindowN ){
+          numNeighbours[idigit2]++;
+        }	   
+      }
+    }
+  }
+
+  // Select by number of neighbours
+  // ==============================
+  for(int idigit=0; idigit<int(DigitList->size()); idigit++){
+    RecoDigit* fdigit = (RecoDigit*)(DigitList->at(idigit));
+    //std::cout << "numNeighbours[" << idigit << "] = " << numNeighbours[idigit] << std::endl;
+    if( numNeighbours[idigit]>=fPmtMinNeighbourDigits ){
+      fSelectByNeighbours->push_back(fdigit);
+    }
+  }
+
+  // delete array of neighbours
+  // ==========================
+  delete [] numNeighbours;
+
+  // return vector of Selected digits
+  // ================================
+  if(verbosity>v_message) std::cout << name << "  Select by neighbours: " << fSelectByNeighbours->size() << std::endl;
+  
+  
+  return fSelectByNeighbours;
+}
+
+std::vector<RecoDigit*>* ClusterSearcher::SelectByClusters(std::vector<RecoDigit*>* DigitList)
+{
+	std::string name = "ClusterSearcher::SelectByClusters() ";
+  // clear vector of Selected digits
+  // ===============================
+  fSelectByClusters->clear();
+  //fRecoClusters->clear();
+
+  // run clustering algorithm
+  // ========================
+  std::vector<RecoCluster>* ClusterList = RecoClusters(DigitList);
+
+  for(int icluster=0; icluster<int(ClusterList->size()); icluster++ ){
+    RecoCluster Cluster = (ClusterList->at(icluster));
+    fRecoClusters->push_back(Cluster);    
+
+    for(int idigit=0; idigit<Cluster.GetNDigits(); idigit++ ){
+      RecoDigit* Digit = new RecoDigit;
+      *Digit=(Cluster.GetDigit(idigit));
+      fSelectByClusters->push_back(Digit);
+    }
+  }
+  
+  // return vector of Selected digits
+  // ================================
+  if(verbosity>v_message) std::cout << name <<"  Select by clusters: " << fSelectByClusters->size() << std::endl;
+  
+  
+  return fSelectByClusters;
+}
+
+std::vector<RecoCluster>* ClusterSearcher::RecoClusters(std::vector<RecoDigit*>* DigitList)
+{  
+
+  // delete cluster digits
+  // =====================
+  for(int i=0; i<int(vClusterDigitList.size()); i++ ){
+    delete (RecoClusterDigit*)(vClusterDigitList.at(i));
+  }
+  vClusterDigitList.clear();
+
+  // clear vector clusters
+  // =====================
+  fClusterList->clear();
+
+  //Digit Selection
+  SelectByPulseHeight(DigitList);
+  SelectByNeighbours(fSelectByPulseHeight);
+
+  // make cluster digits
+  // ===================
+  for(int idigit=0; idigit<int(fSelectByNeighbours->size()); idigit++ ){
+    RecoDigit* recoDigit = (RecoDigit*)(fSelectByNeighbours->at(idigit));
+    RecoClusterDigit* clusterDigit = new RecoClusterDigit(recoDigit);
+    vClusterDigitList.push_back(clusterDigit);
+  }
+
+  // run clustering algorithm
+  // ========================
+  for(int idigit1=0; idigit1<int(vClusterDigitList.size()); idigit1++){
+  	RecoClusterDigit* fdigit1 = (RecoClusterDigit*)(vClusterDigitList.at(idigit1));
+  	int digit1Type = fdigit1->GetDigitType();
+    for(int idigit2=idigit1+1; idigit2<int(vClusterDigitList.size()); idigit2++ ){
+
+      RecoClusterDigit* fdigit2 = (RecoClusterDigit*)(vClusterDigitList.at(idigit2));
+      int digit2Type = fdigit2->GetDigitType();
+
+      double dx = fdigit1->GetX() - fdigit2->GetX();
+      double dy = fdigit1->GetY() - fdigit2->GetY();
+      double dz = fdigit1->GetZ() - fdigit2->GetZ();
+      double dt = fdigit1->GetTime() - fdigit2->GetTime();
+      double drsq = dx*dx + dy*dy + dz*dz;
+      if(digit1Type == RecoDigit::PMT8inch && digit2Type == RecoDigit::PMT8inch) {
+        if( drsq>0.0
+         && drsq<fPmtClusterRadius*fPmtClusterRadius
+         && fabs(dt)<fPmtTimeWindowC ){
+          fdigit1->AddClusterDigit(fdigit2);
+          fdigit2->AddClusterDigit(fdigit1);
+        }	
+        if (drsq == 0 && fabs(dt) < fPmtTimeWindowC) {
+            fdigit2->SetClustered();
+        }
+      }
+      
+      if(digit1Type == RecoDigit::lappd_v0 && digit2Type == RecoDigit::lappd_v0) {
+        if( drsq>0.0
+         && drsq<fLappdClusterRadius*fLappdClusterRadius
+         && fabs(dt)<fLappdTimeWindowC ){
+          fdigit1->AddClusterDigit(fdigit2);
+          fdigit2->AddClusterDigit(fdigit1);
+        }	
+      }
+      
+      if(digit1Type == RecoDigit::PMT8inch && digit2Type == RecoDigit::lappd_v0) {
+        if( drsq>0.0
+         && drsq<fPmtClusterRadius*fPmtClusterRadius
+         && fabs(dt)<fPmtTimeWindowC ){
+          fdigit1->AddClusterDigit(fdigit2);
+        }	
+        if( drsq>0.0
+         && drsq<fLappdClusterRadius*fLappdClusterRadius
+         && fabs(dt)<fLappdTimeWindowC ){
+          fdigit2->AddClusterDigit(fdigit1);
+        }
+      }
+      
+      if(digit1Type == RecoDigit::lappd_v0 && digit2Type == RecoDigit::PMT8inch) {
+        if( drsq>0.0
+         && drsq<fLappdClusterRadius*fLappdClusterRadius
+         && fabs(dt)<fLappdTimeWindowC ){
+          fdigit1->AddClusterDigit(fdigit2);
+        }	
+        if( drsq>0.0
+         && drsq<fPmtClusterRadius*fPmtClusterRadius
+         && fabs(dt)<fPmtTimeWindowC ){
+          fdigit2->AddClusterDigit(fdigit1);
+        }	
+      }
+      
+    }
+  }
+  
+  // collect up clusters
+  // ===================
+  Bool_t carryon = 0;
+  for(int idigit=0; idigit<int(vClusterDigitList.size()); idigit++ ){
+    RecoClusterDigit* fdigit = (RecoClusterDigit*)(vClusterDigitList.at(idigit));
+
+    if( fdigit->IsClustered()==0 && fdigit->GetNClusterDigits()>0 ){
+        Log("ClusterSearcher: valid fdigit index: "+to_string(idigit),v_debug,verbosity);
+        Log("with time: "+to_string(fdigit->GetTime()),v_debug,verbosity);
+        Log("And location: "+to_string(fdigit->GetX())+", "+to_string(fdigit->GetY())+", "+to_string(fdigit->GetZ()),v_debug,verbosity);
+      vClusterDigitCollection.clear();
+      vClusterDigitCollection.push_back(fdigit);
+      fdigit->SetClustered();
+
+      carryon = 1;
+      while( carryon ){
+        carryon = 0;
+        for(int jdigit=0; jdigit<int(vClusterDigitCollection.size()); jdigit++ ){
+	  //std::cout <<"jdigit = "<<jdigit<<", vClusterDigitCollection.size() = "<<vClusterDigitCollection.size()<<std::endl;
+          RecoClusterDigit* cdigit = (RecoClusterDigit*)(vClusterDigitCollection.at(jdigit));
+          TString digitType = cdigit->GetDigitType();
+	        double nDigits = cdigit->GetNClusterDigits();
+	        vNdigitsCluster.push_back(nDigits);	 
+	               
+	        if (digitType==RecoDigit::PMT8inch && cdigit->GetNClusterDigits() >= fPmtMinHitsPerCluster) {
+                
+             if( cdigit->IsAllClustered()==0 ){
+               for( int kdigit=0; kdigit<cdigit->GetNClusterDigits(); kdigit++ ){
+                 RecoClusterDigit* cdigitnew = (RecoClusterDigit*)(cdigit->GetClusterDigit(kdigit));
+                 if( cdigitnew->IsClustered()==0 ){
+                   vClusterDigitCollection.push_back(cdigitnew);
+                   cdigitnew->SetClustered();
+                   carryon = 1;
+                 }
+               }
+             }           
+	        }//end if min of # of hits per cluster
+	        
+	        if (digitType==RecoDigit::lappd_v0 && cdigit->GetNClusterDigits() >= fLappdMinHitsPerCluster) {
+             if( cdigit->IsAllClustered()==0 ){
+               for( int kdigit=0; kdigit<cdigit->GetNClusterDigits(); kdigit++ ){
+                 RecoClusterDigit* cdigitnew = (RecoClusterDigit*)(cdigit->GetClusterDigit(kdigit));
+                 if( cdigitnew->IsClustered()==0 ){
+                   vClusterDigitCollection.push_back(cdigitnew);
+                   cdigitnew->SetClustered();
+                   carryon = 1;
+                 }
+               }
+             }           
+	        }//end if min of # of hits per cluster
+	        
+        }
+      } 
+
+	//std::cout <<"vClusterDigitCollection.size() == "<<vClusterDigitCollection.size()<<std::endl;
+      if( (int)vClusterDigitCollection.size()>=fMinClusterDigits ){
+        RecoCluster cluster;
+        cluster.SetClusterMode(fClusterMode);
+        
+        Log("Adding Digits",v_debug,verbosity);
+        vector<RecoDigit>* ClusteredDigits=new vector<RecoDigit>;
+
+        for(int jdigit=0; jdigit<int(vClusterDigitCollection.size()); jdigit++ ){
+            cout<<"Check 1\n";
+          RecoClusterDigit* cdigit = (RecoClusterDigit*)(vClusterDigitCollection.at(jdigit));
+          cout<<"check 2\n";
+          RecoDigit* arecodigit=cdigit->GetRecoDigit();
+          arecodigit->AddCluster(fClusterMode);
+          RecoDigit brecodigit(arecodigit);
+          ClusteredDigits->push_back(brecodigit);
+          cout<<" check 3\n";
+          //cluster.AddDigit(*arecodigit);        
+        }
+        cluster.SetDigits(ClusteredDigits);
+        ClusteredDigits = nullptr;
+        Log("Cluster has " + to_string(cluster.GetNDigits()) + " digits", v_debug, verbosity);
+        Log("ready to caluclate parameters.",v_debug,verbosity);    
+        cluster.CalcParameters();
+        Log("ClusterSearcher: Clusters made: "+to_string(fClusterList->size()),v_debug,verbosity);
+        fClusterList->push_back(cluster);
+        
+      }
+    }
+  }
+  
+
+  std::cout <<"Number of clusters = "<<fClusterList->size()<<std::endl;
+  // return vector of clusters
+  // =========================
+  return fClusterList;
+}
+
+std::vector<RecoDigit*>* ClusterSearcher::SelectByTruthInfo(std::vector<RecoDigit*>* DigitList)
+{
+	std::string name = " ClusterSearcher::SelectByTruthInfo(() ";
+  // clear vector of Selected digits
+  // ===============================
+  fSelectByTruthInfo->clear();
+  
+  //===================== true info.
+  double x0 = fTrueVertex->GetPosition().X();
+  double y0 = fTrueVertex->GetPosition().Y();
+  double z0 = fTrueVertex->GetPosition().Z();
+  double dirx = fTrueVertex->GetDirection().X();
+  double diry = fTrueVertex->GetDirection().Y();
+  double dirz = fTrueVertex->GetDirection().Z(); 
+  double t0 = 0.0;
+  	
+  // Select by truth
+  // ======================
+  for(int idigit=0; idigit<int(DigitList->size()); idigit++ ){
+    RecoDigit* recoDigit = (RecoDigit*)DigitList->at(idigit);
+    double x = recoDigit->GetPosition().X();
+    double y = recoDigit->GetPosition().Y();
+    double z = recoDigit->GetPosition().Z();
+    double dx = x-x0;
+    double dy = y-y0;
+    double dz = z-z0;
+    double ds = sqrt(dx*dx + dy*dy + dz*dz);
+    double px = dx/ds;
+    double py = dy/ds;
+    double pz = dz/ds;
+    double cosphi = px*dirx + py*diry + pz*dirz;
+    double phideg = acos(cosphi)/3.14*180;
+    if(phideg>38) fSelectByTruthInfo->push_back(recoDigit);   
+  }
+
+  // return vector of Selected digits
+  // ================================
+  if(verbosity>v_message) std::cout << name << "  Select by opening angle: " << fSelectByTruthInfo->size() << std::endl;
+  
+  return fSelectByTruthInfo;
+}
+
diff --git a/UserTools/ClusterSearcher/ClusterSearcher.h b/UserTools/ClusterSearcher/ClusterSearcher.h
new file mode 100644
index 000000000..781e32192
--- /dev/null
+++ b/UserTools/ClusterSearcher/ClusterSearcher.h
@@ -0,0 +1,153 @@
+#ifndef CLUSTERSEARCHER_H
+#define CLUSTERSEARCHER_H
+
+#include <string>
+#include <iostream>
+#include <vector>
+
+#include "Tool.h"
+#include "RecoCluster.h"
+#include "RecoClusterDigit.h"
+#include "TString.h"
+
+class ClusterSearcher: public Tool {
+
+ public:
+
+  ClusterSearcher();
+  ~ClusterSearcher();
+  bool Initialise(std::string configfile,DataModel &data);
+  bool Execute();
+  bool Finalise();
+  
+  typedef enum EFilterConfig {
+    kNone  = 0,
+    kPulseHeight = 1,
+    kPulseHeightAndNeighbours = 2,
+    kPulseHeightAndClusters = 3, 
+    kPulseHeightAndTruthInfo = 4
+  } FilterConfig_t;
+
+  static ClusterSearcher* Instance();
+
+  //static void Config(int config);
+  //static void ClusterType(int ctype);
+  static void PmtMinPulseHeight(double min);
+  static void PmtNeighbourRadius(double radius);
+  static void PmtNeighbourDigits(int digits);
+  static void PmtClusterRadius(double radius);
+  static void MinClusterDigits(int digits);
+  static void PmtTimeWindowN(double windowN);
+  static void PmtTimeWindowC(double windowC);
+
+  void PrintParameters();
+
+  void SetClusterMode(int cmode)		{fClusterMode = cmode;}
+  void SetPmtMinPulseHeight(double min)       { fPmtMinPulseHeight = min; }
+  void SetPmtNeighbourRadius(double radius)   { fPmtNeighbourRadius = radius; }
+  void SetPmtNeighbourDigits(int digits)      { fPmtMinNeighbourDigits = digits; }
+  void SetPmtClusterRadius(double radius)     { fPmtClusterRadius = radius; }
+  void SetPmtTimeWindowNeighbours(double windowN)        { fPmtTimeWindowN = windowN; }
+  void SetPmtTimeWindowClusters(double windowC)        { fPmtTimeWindowC = windowC; }
+  
+  void SetLappdMinPulseHeight(double min)       { fLappdMinPulseHeight = min; }
+  void SetLappdNeighbourRadius(double radius)   { fLappdNeighbourRadius = radius; }
+  void SetLappdNeighbourDigits(int digits)      { fLappdMinNeighbourDigits = digits; }
+  void SetLappdClusterRadius(double radius)     { fLappdClusterRadius = radius; }
+  void SetLappdTimeWindowNeighbours(double windowN)        { fLappdTimeWindowN = windowN; }
+  void SetLappdTimeWindowClusters(double windowC)        { fLappdTimeWindowC = windowC; }
+  //void LoadConfigFile(string configfilename);
+  void SetMinClusterDigits(int digits)        { fMinClusterDigits = digits; }
+  
+
+  std::vector<RecoDigit*>* ResetDigits(std::vector<RecoDigit*>* digitlist);
+  std::vector<RecoDigit*>* SelectDigits(std::vector<RecoDigit*>* digitlist);
+  void SelectDigits(std::vector<RecoDigit>* digitlist);
+  std::vector<RecoDigit*>* SelectAll(std::vector<RecoDigit*>* digitlist);
+  std::vector<RecoDigit*>* SelectByPulseHeight(std::vector<RecoDigit*>* digitlist);
+  std::vector<RecoDigit*>* SelectByNeighbours(std::vector<RecoDigit*>* digitlist);
+  std::vector<RecoDigit*>* SelectByClusters(std::vector<RecoDigit*>* digitlist);
+  std::vector<RecoDigit*>* SelectByTruthInfo(std::vector<RecoDigit*>* digitlist); //use truth information. Only for testing the code
+  std::vector<RecoCluster>* RecoClusters(std::vector<RecoDigit*>* digitlist);
+  
+
+
+ private:
+  void Reset();
+  
+  
+  // clustering mode
+  int fClusterMode;
+
+  // cleaning parameters
+  double fPmtMinPulseHeight;
+  double fPmtNeighbourRadius;
+  int    fPmtMinNeighbourDigits;
+  double fPmtClusterRadius;
+  double fPmtTimeWindowN;
+  double fPmtTimeWindowC;
+  double fPmtMinHitsPerCluster;
+  
+  double fLappdMinPulseHeight;
+  double fLappdNeighbourRadius;
+  int    fLappdMinNeighbourDigits;
+  double fLappdClusterRadius;
+
+  double fLappdTimeWindowN;
+  double fLappdTimeWindowC;
+  double fLappdMinHitsPerCluster;
+  
+  int    fMinClusterDigits;
+  bool   fisMC;
+  bool  fFirstRun;      //Likely temp; would be nice to have it automatically set.
+
+  // p.e. conversion parameters
+  std::map<int,unsigned long> pmt_tubeid_to_channelkey;
+  std::map<unsigned long, double> pmt_gains;
+  std::string singlePEgains;
+
+  // Container for parameters
+  std::map<std::string, double>* fClusteringParam = nullptr;
+
+  // internal containers
+  std::vector<Double_t> vNdigitsCluster;  
+  std::vector<RecoClusterDigit*> vClusterDigitList;
+  std::vector<RecoClusterDigit*> vClusterDigitCollection;
+
+  // vectors of selected digitss
+  std::vector<RecoDigit*>* fSelectAll;
+  std::vector<RecoDigit*>* fSelectByPulseHeight;
+  std::vector<RecoDigit*>* fSelectByNeighbours;
+  std::vector<RecoDigit*>* fSelectByClusters;
+  	
+  // for test only
+  std::vector<RecoDigit*>* fSelectByTruthInfo;
+
+  // vectors of clusters
+  std::vector<RecoCluster>* fClusterList;
+ 
+  // vector of clusters (accessible to the CStore)
+  std::vector<RecoCluster>* fRecoClusters = nullptr;  
+  std::vector<RecoCluster>* pre_RecoClusters = nullptr;
+ 	
+  // true vertex
+  RecoVertex* fTrueVertex = 0; 
+  
+  // digit list
+  std::vector<RecoDigit>* fDigitList = 0;    
+  
+  // hit clustering status;
+  bool fIsHitClusteringDone = false;    
+  
+  /// \brief verbosity levels: if 'verbosity' < this level, the message type will be logged.
+  int verbosity=1;
+	int v_error=0;
+	int v_warning=1;
+	int v_message=2;
+	int v_debug=3;
+	std::string logmessage;
+  
+};
+
+
+#endif
diff --git a/UserTools/ClusterSearcher/README.md b/UserTools/ClusterSearcher/README.md
new file mode 100644
index 000000000..7444b12d7
--- /dev/null
+++ b/UserTools/ClusterSearcher/README.md
@@ -0,0 +1,51 @@
+# ClusterSearcher
+
+ClusterSearcher
+Creates list of RecoClusters according to provided time- and neighboring- parameters.  Clusters are tagged with an int ClusterMode variable for use to differentiate different forms of clusters, created from multiple iterations of the tool within the ToolChain, using separate configfiles (e.g. ClusterSearcherConfig, ClusterSearcherConfig2, etc., using Clustermode 1, 2, etc. respectively)
+
+## Data
+
+Describe any data formats ClusterSearch creates, destroys, changes, or analyzes. E.G.
+
+**pmt_tubeid_to_chanelkey** from CStore
+* Identifies individual PMTs in order to find effective charge in p.e.
+
+**ClusteringParameters** from RecoEvent
+* Compares execution parameters to other iterations of this tool within the toolchain
+
+**TrueVertex** from RecoEvent
+* True vertex position, used for 
+
+**RecoClusters** Into and from RecoEvent
+* List of RecoClusters found based on the configuration parameters provided
+* Any RecoClusters created by previous iteration(s) of ClusterSearcher in the toolchain, so that new clusters can be added to the singular list
+
+
+
+
+## Configuration
+
+
+verbosity                   5
+ClusterMode		    1			#searching track-like clusters	
+Config                      3			#config type: 1 = pulse height cut, 2 = +neighbour cut, 3 = +cluster cut							
+PmtMinPulseHeight           20		#minimum pulse height
+PmtNeighbourRadius          60		#digit neighbouring distance [cm]
+PmtMinNeighbourDigits       2			#minimum neighbour digits
+PmtClusterRadius            60		#digit clustering distance [cm]      
+PmtTimeWindowN              10		#neighbouring time window [ns]    
+PmtTimeWindowC              100		#clustering time window [ns]           
+PmtMinHitsPerCluster        4			#number of hits per cluster								                         
+LappdMinPulseHeight         0			#minimum pulse height                             
+LappdNeighbourRadius        25              #digit neighbouring distance [cm]              
+LappdMinNeighbourDigits     20              #minimum neighbour digits                      
+LappdClusterRadius          25              #digit clustering distance [cm]                                    
+LappdTimeWindowN            1               #neighbouring time window [ns]                 
+LappdTimeWindowC            1               #clustering time window [ns]                   
+LappdMinHitsPerCluster      5			#number of digits per cluster	
+MinClusterDigits            50		#minimum clustered digits	(LAPPD+PMT)		
+MinClusterDigits            4			#minimum clustered digits	(PMT-only)
+SinglePEGains		/path/to/gains/file	#file in which SPE information can be found for PMT charge conversions; only needed on Data	
+FirstRun		    1		#Debug input for differentiating first iteration from later iterations
+
+```
diff --git a/UserTools/DigitBuilder/DigitBuilder.cpp b/UserTools/DigitBuilder/DigitBuilder.cpp
index b3a5a3e4e..3bf94b7f8 100644
--- a/UserTools/DigitBuilder/DigitBuilder.cpp
+++ b/UserTools/DigitBuilder/DigitBuilder.cpp
@@ -18,7 +18,6 @@ DigitBuilder::~DigitBuilder() {
 bool DigitBuilder::Initialise(std::string configfile, DataModel &data){
 
   /////////////////// Usefull header ///////////////////////
-  if(verbosity) cout<<"Initializing Tool DigitBuilder"<<endl;
   if(configfile!="")  m_variables.Initialise(configfile); //loading config file
   //m_variables.Print();
 
@@ -29,7 +28,7 @@ bool DigitBuilder::Initialise(std::string configfile, DataModel &data){
   fParametricModel = 0;
   fIsMC = 1;
   fDigitChargeThr = 10;
-
+  
 
   /// Get the Tool configuration variables
   m_variables.Get("verbosity",verbosity);
@@ -43,9 +42,15 @@ bool DigitBuilder::Initialise(std::string configfile, DataModel &data){
   m_variables.Get("DigitChargeThr",fDigitChargeThr);
   m_variables.Get("ChankeyToPMTIDMap",path_chankeymap);
   m_variables.Get("SinglePEGains",singlePEgains);
+  m_variables.Get("StripHit", striphit);
+  m_variables.Get("MCPMTSmear", MCPMTResolution);
+  m_variables.Get("end_of_window_time_cut", end_of_window_time_cut);
+  m_variables.Get("AcqTimeWindow", AcqTimeWindow);
+  m_variables.Get("CollectHits",fCollectHits);
 
   /// Construct the other objects we'll be setting at event level,
   fDigitList = new std::vector<RecoDigit>;
+  fHitLAPPDs = new std::vector<int>;
 
   // Make the RecoDigit Store if it doesn't exist
   int recoeventexists = m_data->Stores.count("RecoEvent");
@@ -57,6 +62,7 @@ bool DigitBuilder::Initialise(std::string configfile, DataModel &data){
     Log("DigitBuilder Tool: Error retrieving Geometry from ANNIEEvent!",v_error,verbosity); 
     return false; 
   }
+  Log("Strip Hit Mode:" + to_string(striphit),v_debug,verbosity);
   
   // Some hard-coded values of old WCSim LAPPDIDs are in this Tool
   // I would recommend moving away from the use of WCSim IDs if possible as they are liable to change
@@ -66,6 +72,10 @@ bool DigitBuilder::Initialise(std::string configfile, DataModel &data){
     m_data->CStore.Get("channelkey_to_pmtid",channelkey_to_pmtid);
   } else {
     ifstream file_pmtid(path_chankeymap.c_str());
+    if (!file_pmtid) {
+        Log("DigitBuilder Tool: Did not find chankeymap file",v_error,verbosity);
+        return false;
+    }
     while (!file_pmtid.eof()){
       unsigned long chankey;
       int pmtid;
@@ -73,19 +83,39 @@ bool DigitBuilder::Initialise(std::string configfile, DataModel &data){
       channelkey_to_pmtid.emplace(chankey,pmtid);
       pmtid_to_channelkey.emplace(pmtid,chankey);
       if (file_pmtid.eof()) break;
+      Log("DigitBuilder Tool: still gathering chankeys",v_debug,verbosity);
     }
-    
+    Log("DigitBuilder Tool: chankeys done",v_debug,verbosity);
     file_pmtid.close();
     m_data->CStore.Set("pmt_tubeid_to_channelkey",pmtid_to_channelkey);
 
     ifstream file_singlepe(singlePEgains.c_str());
+    if (!file_singlepe) {
+        Log("DigitBuilder Tool: Did not find SinglePEgains file",v_error,verbosity);
+        return false;
+    }
     unsigned long temp_chankey;
     double temp_gain;
-    while (!file_singlepe.eof()){
-      file_singlepe >> temp_chankey >> temp_gain;
-      if (file_singlepe.eof()) break;
-      pmt_gains.emplace(temp_chankey,temp_gain);
+
+    std::string line;
+    if (file_singlepe.is_open()) {
+        //Loop over lines, collect all detector data (should only be one line here)
+        while (getline(file_singlepe, line)) {
+            if (verbosity > 3) std::cout << line << std::endl; //has our stuff;
+            if (line.empty() || line[0] == '#') continue;
+            std::vector<std::string> DataEntries;
+            boost::split(DataEntries, line, boost::is_any_of(","), boost::token_compress_on);
+            int channelkey = -9999;
+            double SPECharge = -9999.;
+            channelkey = std::stoi(DataEntries.at(0));
+            SPECharge = std::stod(DataEntries.at(1));
+            pmt_gains.emplace(channelkey, SPECharge);
+        }
     }
+
+
+    
+    Log("DigitBuilder Tool: SPE gains done",v_debug,verbosity);
     file_singlepe.close();
 
   }
@@ -93,10 +123,10 @@ bool DigitBuilder::Initialise(std::string configfile, DataModel &data){
 
   //Read the LAPPDID file, if given
   if(fLAPPDIDFile!="none"){
-    if(verbosity>2) std::cout << "Loading digits from LAPPD IDs in file " << fLAPPDIDFile << std::endl;
+    Log("Loading digits from LAPPD IDs in file " + fLAPPDIDFile,v_debug,verbosity);
     this->ReadLAPPDIDFile();
   } else {
-    if(verbosity>2) std::cout << "Loading digits from all LAPPDs" << std::endl;
+    Log("Loading digits from all LAPPDs",v_debug,verbosity);
   }
   return true;
 }
@@ -135,16 +165,11 @@ bool DigitBuilder::Execute(){
       return false;
     }
   } else {
-    auto get_clusters =  m_data->CStore.Get("ClusterMap",m_all_clusters);
-    if (!get_clusters){
-      Log("DigitBuilder Tool: ERROR retrieving clustered hits (ClusterMap) in Data mode!",v_error,verbosity);
-      return false;
-    }
-    auto get_clusters_chankey = m_data->CStore.Get("ClusterMapDetkey",m_all_clusters_detkey);
-    if (!get_clusters_chankey){
-      Log("DigitBuilder Tool: ERROR retrieving clustered chankeys (ClusterMapDetkey) in Data mode!",v_error,verbosity);
-      return false;
-    }
+      auto get_dhits = m_data->Stores.at("ANNIEEvent")->Get("Hits",Hits);
+      if (!get_dhits) {
+          Log("DigitBuilder Tool: ERROR retrieving hits in Data mode!",v_error,verbosity);
+          return false;
+      }
   }
 
   /// Build RecoDigit
@@ -160,8 +185,7 @@ bool DigitBuilder::Execute(){
 }
 
 bool DigitBuilder::Finalise(){
-  //delete fDigitList; fDigitList = 0;		//Don't delete pointer to fDigitList, will be deleted by the BoostStore!
-  if(verbosity>0) cout<<"DigitBuilder exitting"<<endl;
+  Log("DigitBuilder exitting",v_message,verbosity);
   return true;
 }
 
@@ -181,7 +205,7 @@ bool DigitBuilder::BuildMCRecoDigit() {
     return true;
   }
   else {
-    cout<<"Wrong PhotoDetector Configuration! Allowed configurations: PMT_only, LAPPD_only, All"<<endl;
+    Log("Wrong PhotoDetector Configuration! Allowed configurations: PMT_only, LAPPD_only, All",v_error);
     return false;
   }
 	
@@ -202,7 +226,7 @@ bool DigitBuilder::BuildDataRecoDigit() {
     return false;
   }
   else {
-    cout<<"Wrong PhotoDetector Configuration! Allowed configurations: PMT_only, LAPPD_only, All"<<endl;
+    Log("Wrong PhotoDetector Configuration! Allowed configurations: PMT_only, LAPPD_only, All",v_error,verbosity);
     return false;
   }
 
@@ -213,13 +237,22 @@ bool DigitBuilder::BuildMCPMTRecoDigit() {
   Log("DigitBuilder Tool: Build PMT reconstructed digits (MC)",v_message,verbosity);
   /// now move to digit retrieval
   int region = -999;
-  double calT;
+  double calT = 0;
   double calQ = 0.;
+  double maxT = -999;
   int digitType = -999;
   Detector* det=nullptr;
   Position  pos_sim, pos_reco;
   /// MCHits is a std::map<unsigned long,std::vector<Hit>>
+
+  
   if(fMCPMTHits){
+      if (fCollectHits) {
+          Log("DigitBuilder Tool: Collecting Hits to make Digits",v_message,verbosity);
+          CollectMCPMTHits();
+          Log("DigitBuilder Tool: # of digits in this event: "+to_string(fDigitList->size()),v_debug,verbosity);
+      }
+      else {
     Log("DigitBuilder Tool: Num PMT Digits = "+to_string(fMCPMTHits->size()),v_message, verbosity);
     /// iterate over the map of sensors with a measurement
     for(std::pair<unsigned long,std::vector<MCHit>>&& apair : *fMCPMTHits){
@@ -243,69 +276,86 @@ bool DigitBuilder::BuildMCPMTRecoDigit() {
       if(det->GetDetectorElement()=="Tank"){
         std::vector<MCHit>& hits = apair.second;
         if(fParametricModel){
-          if(verbosity>2) std::cout << "Using parametric model to build PMT hits" << std::endl;
+          Log("Using parametric model to build PMT hits",v_debug,verbosity);
           //We'll get all hit info and then define a time/charge for each digit
           std::vector<double> hitTimes;
           std::vector<double> hitCharges;
+          //std::vector<int> hitIDs;
           for(MCHit& ahit : hits){
-            if(verbosity>3){
-              std::cout << "This HIT'S TIME AND CHARGE: " << ahit.GetTime() <<
-                  "," << ahit.GetCharge() << std::endl;
-            }
+              Log("This HIT'S TIME AND CHARGE: " + to_string(ahit.GetTime()) + ", " + to_string(ahit.GetCharge()),v_debug,verbosity);
             double hitTime = ahit.GetTime()*1.0;
-          	if(hitTime>-10 && hitTime<40) {
+          	if(hitTime>-10 && hitTime<AcqTimeWindow/1000) {
 			  hitTimes.push_back(ahit.GetTime()*1.0); 
               hitCharges.push_back(ahit.GetCharge());
+              //hitIDs.push_back(ahit.GetHitID());
             }
           }
           // Do median and sum
           std::sort(hitTimes.begin(), hitTimes.end());
           size_t timesize = hitTimes.size();
           if (timesize == 0) continue;
-          if (timesize % 2 == 0){
-            calT = (hitTimes.at(timesize/2 - 1) + hitTimes.at(timesize/2))/2;
-          } else {
-            calT = hitTimes.at(timesize/2);
+          if (fParametricModel == 1) {                //Mean hit time
+              if (timesize % 2 == 0) {
+                  calT = (hitTimes.at(timesize / 2 - 1) + hitTimes.at(timesize / 2)) / 2;
+              }
+              else {
+                  calT = hitTimes.at(timesize / 2);
+              }
           }
+          else if (fParametricModel == 2) {         //first hit time
+              calT = hitTimes.at(0);
+          }
+          else if (fParametricModel == 3) {          //Average of first 20% of hit times
+              for (int i = 0; i < timesize / 5; i++) {
+                  calT += hitTimes.at(i);
+              }
+              if ((int)(timesize / 5) > 0)calT = calT / ((int)(timesize / 5));
+          }
+          else if (fParametricModel == 4) {         //Average hit time
+              for (int i = 0; i < timesize; i++) {
+                  calT += hitTimes.at(i);
+              }
+              calT = calT / timesize;
+          }
+          else if (fParametricModel == 5) {
+              for (int i = 0; i < timesize; i++) {
+                  if (hitCharges.at(i) > maxT) maxT = hitTimes.at(i);
+             }
+              calT = maxT;
+          }
+          if (MCPMTResolution > 0) calT = frand.Gaus(calT, MCPMTResolution);
           calQ = 0.;
           for(std::vector<double>::iterator it = hitCharges.begin(); it != hitCharges.end(); ++it){
             calQ += *it;
           }
-          if (verbosity>4) { 
-            std::cout << "PMT position (X<Y<Z): " << 
-                    to_string(pos_reco.X()) << "," << to_string(pos_reco.Y()) <<
-                    "," << to_string(pos_reco.Z()) << std::endl;
-            std::cout << "PMT Charge,Time: " << to_string(calQ) << "," <<
-                    to_string(calT) << std::endl;
-          }
+            Log("PMT position (X<Y<Z): " + to_string(pos_reco.X()) + "," + to_string(pos_reco.Y()) + "," + to_string(pos_reco.Z()),v_debug,verbosity);
+            Log("PMT Charge,Time: " + to_string(calQ) + "," + to_string(calT),v_debug,verbosity);
+
           if(calQ>fDigitChargeThr) {					//changed to 0 for cross-checks with other tools, change back later!
 				    digitType = RecoDigit::PMT8inch;
 				    RecoDigit recoDigit(region, pos_reco, calT, calQ, digitType, PMTId);
+                    //recoDigit.SetHitIDs(hitIDs);
 				    fDigitList->push_back(recoDigit); 
 				  }
-        } else {
+        }else{
 			    for(MCHit& ahit : hits){
-				  	//if(v_message<verbosity) ahit.Print(); // << VERY verbose
 				  	// get calibrated PMT time (Use the MC time for now)
 				  	calT = ahit.GetTime()*1.0; 
-            calQ = ahit.GetCharge();
-            if (verbosity>4) { 
-              std::cout << "PMT position (X<Y<Z): " << 
-                      to_string(pos_reco.X()) << "," << to_string(pos_reco.Y()) <<
-                      "," << to_string(pos_reco.Z()) << std::endl;
-              std::cout << "PMT Charge,Time: " << to_string(calQ) << "," <<
-                      to_string(calT) << std::endl;
-            }
+            calQ = ahit.GetCharge(); 
+              Log("PMT position (X<Y<Z): " + to_string(pos_reco.X()) + "," + to_string(pos_reco.Y()) + "," + to_string(pos_reco.Z()),v_debug,verbosity);
+              Log("PMT Charge,Time: " + to_string(calQ) + "," + to_string(calT),v_debug,verbosity);
+            calT = frand.Gaus(calT, 1.0);
 				  	digitType = RecoDigit::PMT8inch;
 				  	RecoDigit recoDigit(region, pos_reco, calT, calQ, digitType, PMTId);
-				    //recoDigit.Print();
+                    //recoDigit.SetHitIDs(hitIDs);
 				    fDigitList->push_back(recoDigit); 
           }
-			  }
+        }
       }
 		} // end loop over MCHits
-	} else {
-		cout<<"No MCHits"<<endl;
+    }
+        } else {
+		Log("No MCHits",v_message,verbosity);
 		return false;
 	}
 	return true;
@@ -320,9 +370,14 @@ bool DigitBuilder::BuildMCLAPPDRecoDigit() {
 	int digitType = -999;
 	Detector* det=nullptr;
 	Position  pos_sim, pos_reco;
+    int currentLAPPD = 0;
+    if (fHitLAPPDs->size() > 0) {
+        fHitLAPPDs->clear();
+    }
+
   // repeat for LAPPD hits
 	// MCLAPPDHits is a std::map<unsigned long,std::vector<LAPPDHit>>
-	if(fMCLAPPDHits){
+	if(fMCLAPPDHits && fMCLAPPDHits->size() > 0){
 		Log("DigitBuilder Tool: Num LAPPD Digits = "+to_string(fMCLAPPDHits->size()),v_message,verbosity);
 		// iterate over the map of sensors with a measurement
 		for(std::pair<unsigned long,std::vector<MCLAPPDHit>>&& apair : *fMCLAPPDHits){
@@ -341,46 +396,122 @@ bool DigitBuilder::BuildMCLAPPDRecoDigit() {
       }
       if(!isSelectedLAPPD && fLAPPDId.size()>0) continue;
       if(verbosity>2){
-        std::cout << "Loading in digits for LAPPDID " << LAPPDId << std::endl;
+        Log("Loading in digits for LAPPDID " + to_string(LAPPDId), v_message,verbosity);
+        Log("located: ",v_message,verbosity);
+        det->GetPositionInTank().Print();
+        Log("directed: ",v_message,verbosity);
+        det->GetDetectorDirection().Print();
       }
 
-			if(det->GetDetectorElement()=="LAPPD"){ // redundant, MCLAPPDHits are LAPPD hitss
-				std::vector<MCLAPPDHit>& hits = apair.second;
-				for(MCLAPPDHit& ahit : hits){
-					//if(v_message<verbosity) ahit.Print(); // << VERY verbose
-					// an LAPPDHit has adds (global x-y-z) position, (in-tile x-y) local position
-					// and time psecs
-					// convert the WCSim coordinates to the ANNIEreco coordinates
-					// convert the unit from m to cm
-					pos_reco.SetX(ahit.GetPosition().at(0)*100.+xshift); //cm
-					pos_reco.SetY(ahit.GetPosition().at(1)*100.+yshift); //cm
-					pos_reco.SetZ(ahit.GetPosition().at(2)*100.+zshift); //cm
-					calT = ahit.GetTime();  // 
-					calT = frand.Gaus(calT, 0.1); // time is smeared with 100 ps time resolution. Harded-coded for now.
-					calQ = ahit.GetCharge();
-          if (verbosity>4) { 
-            std::cout << "LAPPD position (X<Y<Z): " << 
-                    to_string(pos_reco.X()) << "," << to_string(pos_reco.Y()) <<
-                    "," << to_string(pos_reco.Z()) << std::endl;
-            std::cout << "LAPPD Charge,Time: " << to_string(calQ) << "," <<
-                    to_string(calT) << std::endl;
+      if (det->GetDetectorElement() == "LAPPD") { // redundant, MCLAPPDHits are LAPPD hitss
+          std::vector<MCLAPPDHit>& hits = apair.second;
+
+          std::vector<double> sumposX;
+          std::vector<double> sumposY;
+          std::vector<double> sumposZ;
+          std::vector<double> sumT;
+          double posX;
+          std::vector<int> nHitsOnStrip;
+          int a;
+          for (int i = 0; i < 28; i++) {
+              sumposX.push_back(0);
+              sumposY.push_back(0);
+              sumposZ.push_back(0);
+              sumT.push_back(0);
+              nHitsOnStrip.push_back(0);
           }
-					// I found the charge is 0 for all the hits. In order to test the code, 
-					// here I just set the charge to 1. We should come back to this later. (Jingbo Wang)
-					calQ = 1.;
-					digitType = RecoDigit::lappd_v0;
-					RecoDigit recoDigit(region, pos_reco, calT, calQ, digitType,LAPPDId);
-					//if(v_message<verbosity) recoDigit.Print();
-				  //make some cuts here. It will be moved to the Hitcleaning tool
-				  if(calT>40 || calT<-10) continue; // cut off delayed hits
-				  fDigitList->push_back(recoDigit);
+				for(MCLAPPDHit& ahit : hits){
+                    if (LAPPDId != currentLAPPD && hits.size() > 3) {
+                        Log("THERE ARE HITS ON LAPPD " + to_string(LAPPDId),v_message,verbosity);
+                        currentLAPPD = LAPPDId;                        
+                        fHitLAPPDs->push_back(currentLAPPD);
+                        Log("fHitLAPPDs size: " + to_string(fHitLAPPDs->size()),v_message,verbosity);
+                    }
+                    if (striphit == 0) {
+                        //if(v_message<verbosity) ahit.Print(); // << VERY verbose
+                        // an LAPPDHit has adds (global x-y-z) position, (in-tile x-y) local position
+                        // and time psecs
+                        // convert the WCSim coordinates to the ANNIEreco coordinates
+                        // convert the unit from m to cm
+                        pos_reco.SetX(ahit.GetPosition().at(0) * 100. + xshift); //cm
+                        pos_reco.SetY(ahit.GetPosition().at(1) * 100. + yshift); //cm
+                        pos_reco.SetZ(ahit.GetPosition().at(2) * 100. + zshift); //cm
+                        calT = ahit.GetTime();  // 
+                        calT = frand.Gaus(calT, 0.1); // time is smeared with 100 ps time resolution. Harded-coded for now.
+                        calQ = ahit.GetCharge();
+                            
+                        Log("LAPPD position (X<Y<Z): " + to_string(pos_reco.X()) + "," + to_string(pos_reco.Y()) + "," + to_string(pos_reco.Z()),v_debug,verbosity);
+                            
+                        Log("LAPPD Charge,Time: " + to_string(calQ) + "," + to_string(calT),v_debug,verbosity);
+                        // I found the charge is 0 for all the hits. In order to test the code, 
+                        // here I just set the charge to 1. We should come back to this later. (Jingbo Wang)
+                        calQ = 1.;
+                        digitType = RecoDigit::lappd_v0;
+                        RecoDigit recoDigit(region, pos_reco, calT, calQ, digitType, LAPPDId);
+                        //if(v_message<verbosity) recoDigit.Print();
+                      //make some cuts here. It will be moved to the Hitcleaning tool
+                        if (calT > 40 || calT < -10) continue; // cut off delayed hits
+                        //recoDigit.SetHitIDs(hitIDs);
+                        fDigitList->push_back(recoDigit);
+
+                    }
+                    else {
+                        posX = (ahit.GetPosition().at(0) * 100 + xshift) - (det->GetPositionInTank().X() * 100 + xshift);
+                        double dirX = det->GetDetectorDirection().X();
+                        double AngX = std::asin(dirX);
+                        int strip = (int)((posX / (20.* std::sin(AngX) / 28.)) + 14 - 1);
+                        if (strip > 27 || strip < 0) {
+                            Log("warning: LAPPD hit found off LAPPD: ", v_warning, verbosity);
+                            Log("hit found on strip " + to_string(strip),v_warning,verbosity);
+                            Log("LAPPD truehit X,Y,Z: " + to_string(ahit.GetPosition().at(0)) + ", " + to_string(ahit.GetPosition().at(1)) + ", " + to_string(ahit.GetPosition().at(2)),v_warning,verbosity);
+                            Log("continuing loop without this hit.",v_warning,verbosity);
+                            continue;
+                        }
+                        Log("hit found on strip " + to_string(strip),v_debug,verbosity);
+                        sumposX.at(strip) += ahit.GetPosition().at(0) * 100 + xshift;
+                        sumposY.at(strip) += ahit.GetPosition().at(1) * 100 + yshift;
+                        sumposZ.at(strip) += ahit.GetPosition().at(2) * 100 + zshift;
+                        Log("LAPPD truehit X,Y,Z: " + to_string(ahit.GetPosition().at(0)) + ", " + to_string(ahit.GetPosition().at(1)) + ", " + to_string(ahit.GetPosition().at(2)),v_debug,verbosity);
+                        
+                        if ((striphit == 3 && nHitsOnStrip.at(strip) < 3) || (striphit==2 && nHitsOnStrip.at(strip) == 0) || striphit == 1) 
+                            sumT.at(strip) += frand.Gaus(calT, 0.1); // time is smeared with 100 ps time resolution. Harded-coded for now.
+                        nHitsOnStrip.at(strip) += 1;
+
+                    }
 				}
+                if (striphit != 0) {
+                    for (int strip = 0; strip < 28; strip++) {
+                        if (nHitsOnStrip.at(strip) != 0) {
+                            pos_reco.SetX(sumposX.at(strip) / nHitsOnStrip.at(strip));
+                            pos_reco.SetY(sumposY.at(strip) / nHitsOnStrip.at(strip));
+                            pos_reco.SetZ(sumposZ.at(strip) / nHitsOnStrip.at(strip));
+                            if(striphit == 1) calT = sumT.at(strip) / nHitsOnStrip.at(strip);
+                            if (striphit == 3) calT = sumT.at(strip) / 3;
+                            calQ = nHitsOnStrip.at(strip); //set to the number of hits for now.
+                            digitType = RecoDigit::lappd_v0;
+
+                            Log("LAPPD ID: " + to_string(LAPPDId),v_message,verbosity);
+                            Log("LAPPD strip-hit position (X<Y<Z): " + to_string(pos_reco.X()) + "," + to_string(pos_reco.Y()) + "," + to_string(pos_reco.Z()),v_debug,verbosity);
+                            Log("LAPPD strip-hit Charge,Time: " + to_string(calQ) + "," + to_string(calT),v_debug,verbosity);
+                            RecoDigit recoDigit(region, pos_reco, calT, calQ, digitType, LAPPDId);
+                            Log("recoDigit created for striphit",v_message,verbosity);
+                            //recoDigit.SetHitIDs(hitIDs);
+                            fDigitList->push_back(recoDigit);
+                        }
+                    }
+                }
+                sumposX.clear();
+                sumposY.clear();
+                sumposZ.clear();
+                sumT.clear();
+                nHitsOnStrip.clear();
 			}
 		} // end loop over MCLAPPDHits
 	} else {
-		cout<<"No MCLAPPDHits"<<endl;
+		Log("No MCLAPPDHits",v_warning,verbosity);
 		return false;
 	}
+    
 	return true;
 }
 
@@ -388,161 +519,43 @@ bool DigitBuilder::BuildDataPMTRecoDigit(){
 
 	Log("DigitBuilder Tool: Build PMT reconstructed digits (data)",v_message,verbosity);
 	/// now move to digit retrieval
-	int region = -999;
-	double calT;
-	double calQ = 0.;
-	int digitType = -999;
-	Detector* det=nullptr;
-	Position  pos_sim, pos_reco;
-	/// m_all_clusters is a std::map<double,std::vector<Hit>>
-        
-	if (m_all_clusters && m_all_clusters_detkey){
-          int clustersize = m_all_clusters->size();
-          std::cout <<"Clustersize of m_all_clusters: "<<clustersize<<std::endl;
-          bool clusters_available = false;
-          bool muon_available = false;
-          if (clustersize != 0) clusters_available = true;
-          if (clusters_available){
-	  //determine the main cluster (max charge and in [0 ... 2000ns] time window)
-	  double max_cluster = 0;
-          double max_charge = 0;
-          for(std::pair<double,std::vector<Hit>>&& apair : *m_all_clusters){
-            std::vector<Hit>&Hits = apair.second;
-            double time = 0;
-            int hits=0;
-            double charge = 0;
-            for (unsigned int i_hit = 0; i_hit < Hits.size(); i_hit++){
-              hits++;
-              time+=Hits.at(i_hit).GetTime();
-              charge+=Hits.at(i_hit).GetCharge();
-          }
-          if (hits>0) {
-            time/=hits;
-          }
-          if (time > 2000.) continue;	//not a beam muon if not in primary window
-	  if (charge > max_charge) {
-            muon_available = true;
-            max_charge = charge;
-            max_cluster = apair.first;
-          }
-	}
-	if (muon_available){
-	  std::vector<Hit>& Hits = m_all_clusters->at(max_cluster);
-          std::vector<unsigned long> detkeys = m_all_clusters_detkey->at(max_cluster);
-          int hits_pmt = 0;
-
-          std::map<unsigned long,std::vector<double>> hitTimes;
-          std::map<unsigned long,std::vector<double>> hitCharges;
-
-	  Log("DigitBuilder Tool: Num PMT Clustered Digits = "+to_string(Hits.size()),v_message, verbosity);
-	  for (unsigned int i_hit = 0; i_hit < Hits.size(); i_hit++){
-	    Hit ahit = Hits.at(i_hit);
-            unsigned long chankey = detkeys.at(i_hit);
-	    
-
-	    if (hitTimes.find(chankey)!=hitTimes.end()){
-	      hitTimes.at(chankey).push_back(ahit.GetTime());
-	      hitCharges.at(chankey).push_back(ahit.GetCharge());
-            }
-	    else {
-	      std::vector<double> temp_hittimes{ahit.GetTime()};
-	      std::vector<double> temp_hitcharges{ahit.GetCharge()};
-              hitTimes.emplace(chankey,temp_hittimes);
-              hitCharges.emplace(chankey,temp_hitcharges);
-	    }
-          }
+    int region = -999;
+    double calT;
+    double calQ = 0.;
+    double calQ_temp = 0;
+    int digitType = -999;
+    Detector* det = nullptr;
+    Position  pos_sim, pos_reco;
+
+    for (std::pair<unsigned long, std::vector<Hit>>&& apair : *Hits) {
+        unsigned long chankey = apair.first;
+        Detector* thistube = fGeometry->ChannelToDetector(chankey);
+        det = fGeometry->ChannelToDetector(chankey);
+        int detectorkey = thistube->GetDetectorID();
+        int PMTId = channelkey_to_pmtid.at(chankey);
+        if (thistube->GetDetectorElement() == "Tank") {
+            std::vector<Hit>& ThisPMTHits = apair.second;
+            for (Hit& ahit : ThisPMTHits) {
+                pos_reco=det->GetDetectorPosition();
+                calT = ahit.GetTime();
+
+                calQ = ahit.GetCharge();
+          
+                if (pmt_gains.find(chankey) != pmt_gains.end() && pmt_gains.at(chankey) > 0.0) {
+                    calQ_temp = calQ / pmt_gains.at(chankey);
+                }
+                if (calQ_temp > fDigitChargeThr) {
+                    digitType = RecoDigit::PMT8inch;
+                    RecoDigit recoDigit(region, pos_reco, calT, calQ_temp, digitType, PMTId);
 
-          if(fParametricModel){
-            Log("DigitBuilder tool: Use Parametric Model to create digits",v_message,verbosity);
-            // Do median and sum
-            std::map<unsigned long,std::vector<double>>::iterator it, it2;
-            for (it=hitTimes.begin(),it2 = hitCharges.begin(); it != hitTimes.end(), it2 != hitCharges.end(); it++, it2++){
-	      unsigned long chankey = it->first;
-	      std::vector<double> hittimes = it->second;
-	      std::vector<double> hitcharges = it2->second;
-	      det = fGeometry->ChannelToDetector(chankey);
-	      int PMTId = channelkey_to_pmtid.at(chankey);  //PMTID In WCSim
-	      if(det==nullptr){
-	        Log("DigitBuilder Tool: Detector not found! ",v_message,verbosity);
-	        continue;
-	      }
-              // convert the WCSim coordinates to the ANNIEreco coordinates
-	      // convert the unit from m to cm
-	      pos_sim = det->GetDetectorPosition();
-	      pos_sim.UnitToCentimeter();
-	      pos_reco.SetX(pos_sim.X()+xshift);
-	      pos_reco.SetY(pos_sim.Y()+yshift);
-	      pos_reco.SetZ(pos_sim.Z()+zshift);
-
-              std::sort(hittimes.begin(), hittimes.end());
-              size_t timesize = hittimes.size();
-              if (timesize == 0) continue;
-              if (timesize % 2 == 0){
-                calT = (hittimes.at(timesize/2 - 1) + hittimes.at(timesize/2))/2;
-              } else {
-                calT = hittimes.at(timesize/2);
-              }
-              calQ = 0.;
-              for(std::vector<double>::iterator it3 = hitcharges.begin(); it3 != hitcharges.end(); ++it3){
-                calQ += *it3;
-              }
-              if (verbosity>4) { 
-                std::cout << "PMT position (X<Y<Z): " << 
-                to_string(pos_reco.X()) << "," << to_string(pos_reco.Y()) <<
-                "," << to_string(pos_reco.Z()) << std::endl;
-                std::cout << "PMT Charge,Time: " << to_string(calQ) << "," <<
-                to_string(calT) << std::endl;
-              }
-	      double calQ_temp = calQ;
-	      if (pmt_gains.find(chankey) != pmt_gains.end() && pmt_gains.at(chankey) > 0.0){
-                calQ_temp = calQ / pmt_gains.at(chankey);
-              }
-              if(calQ_temp>fDigitChargeThr) {
-                digitType = RecoDigit::PMT8inch;
-	        RecoDigit recoDigit(region, pos_reco, calT, calQ_temp, digitType, PMTId);
-                fDigitList->push_back(recoDigit); 
-	      }
-            }
-          } else {
-            std::map<unsigned long,std::vector<double>>::iterator it, it2;
-            for (it=hitTimes.begin(), it2 = hitCharges.begin(); it != hitTimes.end(), it2 != hitCharges.end(); it++, it2++){
-              unsigned long chankey = it->first;
-              std::vector<double> hittimes = it->second;
-              std::vector<double> hitcharges = it2->second;
-              det = fGeometry->ChannelToDetector(chankey);
-              int PMTId = channelkey_to_pmtid.at(chankey);  //PMTID In WCSim
-              if(det==nullptr){
-                Log("DigitBuilder Tool: Detector not found! ",v_message,verbosity);
-                continue;
-              }
-              pos_sim = det->GetDetectorPosition();
-              pos_sim.UnitToCentimeter();
-              pos_reco.SetX(pos_sim.X()+xshift);
-              pos_reco.SetY(pos_sim.Y()+yshift);
-              pos_reco.SetZ(pos_sim.Z()+zshift);
-                
-              for (int i=0; i< int(hitcharges.size()); i++){
-		calT = hittimes.at(i);
-                calQ = hitcharges.at(i);
-                if (verbosity>4) {
-                  std::cout << "PMT position (X<Y<Z): " <<
-                  to_string(pos_reco.X()) << "," << to_string(pos_reco.Y()) <<
-                  "," << to_string(pos_reco.Z()) << std::endl;
-                  std::cout << "PMT Charge,Time: " << to_string(calQ) << "," <<
-                  to_string(calT) << std::endl;
+                    //recoDigit.SetHitIDs(hitIDs);
+                    fDigitList->push_back(recoDigit);
                 }
-                digitType = RecoDigit::PMT8inch;
-                RecoDigit recoDigit(region, pos_reco, calT, calQ, digitType, PMTId);
-		fDigitList->push_back(recoDigit);
-               }
-             }
-           } 
-         }
-       }
-     } else {
-       cout<<"No Clustered Hits found."<<endl;
-       return false;
-     }
+
+            }
+            
+        }
+    }
      
      return true;
 
@@ -551,11 +564,13 @@ bool DigitBuilder::BuildDataPMTRecoDigit(){
 void DigitBuilder::PushRecoDigits(bool savetodisk) {
 	Log("DigitBuilder Tool: Push reconstructed digits to the RecoEvent store",v_message,verbosity);
 	m_data->Stores.at("RecoEvent")->Set("RecoDigit", fDigitList, savetodisk);  ///> Add digits to RecoEvent
+    m_data->Stores.at("RecoEvent")->Set("HitLAPPDs", fHitLAPPDs, savetodisk);
 }
 
 void DigitBuilder::Reset() {
   // Reset 
   fDigitList->clear();
+  fHitLAPPDs->clear();
 }
 
 void DigitBuilder::ReadLAPPDIDFile() {
@@ -564,7 +579,7 @@ void DigitBuilder::ReadLAPPDIDFile() {
   if (myfile.is_open()){
     while(getline(myfile,line)){
       if(verbosity>0){
-        std::cout << "DigitBuilder tool: Loading hits from LAPPD ID " << line << std::endl;
+        Log("DigitBuilder tool: Loading hits from LAPPD ID " + line,v_message,verbosity);
       }
       int thisID = std::atoi(line.c_str());
       fLAPPDId.push_back(thisID);
@@ -573,3 +588,145 @@ void DigitBuilder::ReadLAPPDIDFile() {
     Log("Unable to open given LAPPD ID File. Using all LAPPDs",v_error,verbosity);
   }
 }
+
+void DigitBuilder::CollectMCPMTHits() {
+    double calT = 0;
+    double calQ = 0;
+    Position pos_reco, pos_sim;
+    std::map<unsigned long, int> PMT_ishit;
+
+
+    for (std::pair<unsigned long, std::vector<MCHit>>&& apair : *fMCPMTHits) {
+        unsigned long chankey = apair.first;
+        Detector* thistube = fGeometry->ChannelToDetector(chankey);
+        int detectorkey = thistube->GetDetectorID();
+        int PMTId = channelkey_to_pmtid.at(chankey);
+        Log("DigitBuilder Tool: collecting hits for PMT "+to_string(detectorkey)+", which corresponds to PMTId "+to_string(PMTId),v_debug,verbosity);
+        if (thistube->GetDetectorElement() == "Tank") {
+            std::vector<MCHit>& hits = apair.second;
+            PMT_ishit[detectorkey] = 1;
+            std::vector<double> hit_times;
+            std::vector<double> temp_times;
+            double first_time=0;
+            double mid_time=0;
+            std::vector<double> datalike_hits;
+            std::vector<double> datalike_hits_charge;
+            std::vector<int> Parents_by_hit;
+            std::vector<int> temp_parents;
+            double max_charge=-999;
+            pos_sim = thistube->GetDetectorPosition();
+            pos_sim.UnitToCentimeter();
+            pos_reco.SetX(pos_sim.X() + xshift);
+            pos_reco.SetY(pos_sim.Y() + yshift);
+            pos_reco.SetZ(pos_sim.Z() + zshift);
+
+
+            for (MCHit& ahit : hits) {
+
+                    hit_times.push_back(ahit.GetTime());
+
+            }
+            
+
+            if (hit_times.size() == 0) {
+                Log("DigitBuilder tool: no hits in window.",v_message,verbosity);
+                return;
+            }
+
+            //Combine multiple MC hits to one pulse
+            std::sort(hit_times.begin(), hit_times.end());
+            for (int i_hit = 0; i_hit < (int)hit_times.size(); i_hit++) {
+                Log("Hits, hit_times, datalike_hits, parents size: "+to_string(hits.size())+" "+to_string(hit_times.size())+" "+to_string(datalike_hits.size())+" " + to_string(hits.at(i_hit).GetParents()->size()), v_debug, verbosity);
+                double hit1 = hit_times.at(i_hit);
+                if(hit1>end_of_window_time_cut * AcqTimeWindow) continue;
+                int j_hit=0;
+                if (datalike_hits.size() == 0) {
+
+                    
+                    datalike_hits.push_back(hit1);
+                    first_time=hit1;
+
+                    datalike_hits_charge.push_back(hits.at(i_hit).GetCharge());
+                    if (hits.at(i_hit).GetParents()->size() == 0) {
+                        Parents_by_hit.push_back(-5);
+                        Log("found hit with no parent at time "+to_string(hits.at(i_hit).GetTime()),v_debug,verbosity);
+                    }
+                    else {
+                        Parents_by_hit.push_back(hits.at(i_hit).GetParents()->at(0));
+                    }
+
+                    temp_times.push_back(hit1);
+                }
+                else {
+                    bool new_pulse = false;
+
+                    for (int j_hit = 0; j_hit < (int)datalike_hits.size(); j_hit++) {
+                        
+                        if (fabs(first_time - hit1) < 10.) {
+                            new_pulse = false;
+                            datalike_hits_charge.at(j_hit) += hits.at(i_hit).GetCharge();
+                            temp_times.push_back(hit1);
+                            
+                        }
+                        else {
+                            new_pulse = true;
+                            temp_times.push_back(hit1);
+                            
+                        }
+
+                    }
+
+                    if (new_pulse) {
+                        first_time=hit1;
+                        
+                        // following the DigitBuilder tool --> take median photon hit time as the hit time of the "pulse"
+                        sort(temp_times.begin(),temp_times.end());
+                        if (temp_times.size() % 2 == 0) {
+                            mid_time = (temp_times.at(temp_times.size() / 2 - 1) + temp_times.at(temp_times.size() / 2)) / 2;
+                        }
+                        else {
+                            mid_time = temp_times.at(temp_times.size() / 2);
+                        }
+
+                        datalike_hits.at(j_hit)=mid_time;  //datalike_hits.push_back(mid_time);      //Only count as a new pulse if it was 10ns away from every other pulse
+                        datalike_hits_charge.push_back(hits.at(i_hit).GetCharge());
+                        j_hit++;
+
+                        datalike_hits.push_back(hit1);
+
+                        if (hits.at(i_hit).GetParents()->size() == 0) Parents_by_hit.push_back(-5);
+                        else {
+                            Parents_by_hit.push_back(hits.at(i_hit).GetParents()->at(0));
+                        Log("Hit parent " + to_string(hits.at(i_hit).GetParents()->at(0)) + " and time: " + to_string(hits.at(i_hit).GetTime()), v_debug, verbosity);
+                        }
+                        temp_times.clear();
+
+                    }
+                }
+            }
+
+
+            //hits.clear();
+            
+            for (int i_hit = 0; i_hit < (int)datalike_hits.size(); i_hit++) {
+
+                calT = datalike_hits.at(i_hit);
+
+                if (MCPMTResolution > 0) calT= frand.Gaus(calT, MCPMTResolution);
+                calQ = datalike_hits_charge.at(i_hit);
+
+                RecoDigit recoDigit(-999 /*region*/, pos_reco, calT, calQ, RecoDigit::PMT8inch, PMTId);
+                temp_parents.push_back(Parents_by_hit.at(i_hit));
+
+                recoDigit.SetParents(temp_parents);
+                fDigitList->push_back(recoDigit);
+                temp_parents.clear();
+
+                Log("PMT position (X<Y<Z): " + to_string(pos_reco.X()) + "," + to_string(pos_reco.Y()) + "," + to_string(pos_reco.Z()), v_debug, verbosity);
+                Log("PMT Charge,Time: " + to_string(calQ) + "," + to_string(calT), v_debug, verbosity);
+                //newMCHits.push_back(MCHit(chankey, datalike_hits.at(i_hit), datalike_hits_charge.at(i_hit), parents));
+            }
+            
+        }
+    }
+}
\ No newline at end of file
diff --git a/UserTools/DigitBuilder/DigitBuilder.h b/UserTools/DigitBuilder/DigitBuilder.h
index ed3986175..18013a722 100644
--- a/UserTools/DigitBuilder/DigitBuilder.h
+++ b/UserTools/DigitBuilder/DigitBuilder.h
@@ -16,6 +16,7 @@
 #include "TTree.h"
 #include "ANNIEGeometry.h"
 #include "Detector.h"
+#include <boost/algorithm/string.hpp>
 
 class DigitBuilder: public Tool {
 
@@ -43,6 +44,7 @@ class DigitBuilder: public Tool {
   /// It adds PMT hits to the RecoDigit list
    bool BuildMCPMTRecoDigit();
    bool BuildDataPMTRecoDigit();   ///> Same as BuildMCPMTRecoDigit, but applicable for data.
+   void CollectMCPMTHits();
  	
   /// \brief Build LAPPD digits
   ///
@@ -72,15 +74,21 @@ class DigitBuilder: public Tool {
   int verbosity=1;
   std::string fInputfile;
   unsigned long fNumEvents;
-  
+  std::vector<int>* fHitLAPPDs;
   std::vector<int> fLAPPDId; ///< selected LAPPDs
   std::string fPhotodetectorConfiguration; ///< "PMTs_Only", "LAPPDs_Only", "All_Detectors"
-  bool fParametricModel;     ///< configures if PMTs hits for each event are accumulated into one hit per PMT
+  int fParametricModel;     ///< configures how PMTs hits for each event are accumulated into one hit per PMT 0: they are not, 1: median hit time, 2: first hit time, 3: average first 20% of hits, 4: average all hits, 5: highest charge hit time
   bool fIsMC;     ///< Configure whether to load from MCHits or Hits in boost store 
   std::string  fLAPPDIDFile="none";
   double fDigitChargeThr;
   std::string path_chankeymap;
   std::string singlePEgains;
+  int striphit; //0 all LAPPD Hits as true; 1 average all times per strip; 2 use first time for each strip
+  double MCPMTResolution = 0;
+
+  bool fCollectHits;
+  int AcqTimeWindow;
+  double end_of_window_time_cut;
 
   Geometry* fGeometry=nullptr;    ///< ANNIE Geometry
   TRandom3 frand;  ///< Random number generator
@@ -104,10 +112,11 @@ class DigitBuilder: public Tool {
   /// Reconstructed information
   std::vector<RecoDigit>* fDigitList;				///< Reconstructed Hits including both LAPPD hits and PMT hits
   std::map<unsigned long,std::vector<MCHit>>* fMCPMTHits=nullptr;             ///< PMT hits
+  std::map<unsigned long, std::vector<Hit>>* Hits = nullptr;
   std::map<unsigned long,std::vector<MCLAPPDHit>>* fMCLAPPDHits=nullptr;   ///< LAPPD hits
   std::map<unsigned long,std::vector<MCHit>>* fTDCData=nullptr;            ///< MRD & veto hits
-  std::map<double,std::vector<Hit>>* m_all_clusters=nullptr;            ///< Clusters, from ClusterFinder tool
-  std::map<double,std::vector<unsigned long>>* m_all_clusters_detkey=nullptr;         ///< Chankeys corresponding to clusters, from ClusterFinder tool
+  std::map<double,std::vector<Hit>>* m_all_clusters=nullptr;            ///< Clusters, from ClusterFinder tool - deprecated and to be removed
+  std::map<double,std::vector<unsigned long>>* m_all_clusters_detkey=nullptr;         ///< Chankeys corresponding to clusters, from ClusterFinder tool - deprecated and to be removed
 
   std::map<unsigned long, double> pmt_gains;
 
diff --git a/UserTools/Factory/Factory.cpp b/UserTools/Factory/Factory.cpp
index 5383362c2..b73507d4d 100644
--- a/UserTools/Factory/Factory.cpp
+++ b/UserTools/Factory/Factory.cpp
@@ -178,6 +178,8 @@ if (tool=="EBPMT") ret=new EBPMT;
 if (tool=="LAPPDLoadStore") ret=new LAPPDLoadStore;
 if (tool=="PMTWaveformSim") ret=new PMTWaveformSim;
 if (tool=="LAPPDWaveformDisplay") ret=new LAPPDWaveformDisplay;
+if (tool=="ClusterSearcher") ret=new ClusterSearcher;
+if (tool=="NeutronCheck") ret=new NeutronCheck;
 if (tool=="PrintADCTraces") ret=new PrintADCTraces;
 return ret;
 }
diff --git a/UserTools/HitCleaner/HitCleaner.cpp b/UserTools/HitCleaner/HitCleaner.cpp
index 3475b4ce2..4adc803d9 100644
--- a/UserTools/HitCleaner/HitCleaner.cpp
+++ b/UserTools/HitCleaner/HitCleaner.cpp
@@ -118,8 +118,8 @@ bool HitCleaner::Initialise(std::string configfile, DataModel &data){
   // only for test 
   fFilterByTruthInfo = new std::vector<RecoDigit*>; 
   // vector of clusters
-  fClusterList = new std::vector<RecoCluster*>;
-  fHitCleaningClusters = new std::vector<RecoCluster*>;
+  fClusterList = new std::vector<RecoCluster>;
+  fHitCleaningClusters = new std::vector<RecoCluster>;
 
   //Set hit cleaner parameters in the RecoEvent store
   m_data->Stores.at("RecoEvent")->Set("HitCleaningParameters", fHitCleaningParam);
@@ -194,6 +194,7 @@ bool HitCleaner::Execute(){
   fIsHitCleaningDone = true;
   m_data->Stores.at("RecoEvent")->Set("HitCleaningDone", fIsHitCleaningDone); 
   m_data->Stores.at("RecoEvent")->Set("HitCleaningClusters", fHitCleaningClusters);
+  CBCheck(digits,FilterDigitList);
 
   delete digits; digits = 0;
   return true;
@@ -523,14 +524,15 @@ std::vector<RecoDigit*>* HitCleaner::FilterByClusters(std::vector<RecoDigit*>* m
 
   // run clustering algorithm
   // ========================
-  std::vector<RecoCluster*>* myClusterList = (std::vector<RecoCluster*>*)(this->RecoClusters(myDigitList));
+  std::vector<RecoCluster>* myClusterList = (std::vector<RecoCluster>*)(this->RecoClusters(myDigitList));
 
   for(int icluster=0; icluster<int(myClusterList->size()); icluster++ ){
-    RecoCluster* myCluster = (RecoCluster*)(myClusterList->at(icluster));
+    RecoCluster myCluster = (myClusterList->at(icluster));
     fHitCleaningClusters->push_back(myCluster);    
 
-    for(int idigit=0; idigit<myCluster->GetNDigits(); idigit++ ){
-      RecoDigit* myDigit = (RecoDigit*)(myCluster->GetDigit(idigit));
+    for(int idigit=0; idigit<myCluster.GetNDigits(); idigit++ ){
+      RecoDigit* myDigit = new RecoDigit;
+      *myDigit=(myCluster.GetDigit(idigit));
       fFilterByClusters->push_back(myDigit);
     }
   }
@@ -543,7 +545,7 @@ std::vector<RecoDigit*>* HitCleaner::FilterByClusters(std::vector<RecoDigit*>* m
   return fFilterByClusters;
 }
 
-std::vector<RecoCluster*>* HitCleaner::RecoClusters(std::vector<RecoDigit*>* myDigitList)
+std::vector<RecoCluster>* HitCleaner::RecoClusters(std::vector<RecoDigit*>* myDigitList)
 {  
 
   // delete cluster digits
@@ -679,13 +681,13 @@ std::vector<RecoCluster*>* HitCleaner::RecoClusters(std::vector<RecoDigit*>* myD
       } 
 	//std::cout <<"vClusterDigitCollection.size() == "<<vClusterDigitCollection.size()<<std::endl;
       if( (int)vClusterDigitCollection.size()>=fMinClusterDigits ){
-        RecoCluster* cluster = new RecoCluster();
+        RecoCluster cluster;
         fClusterList->push_back(cluster);
 
         for(int jdigit=0; jdigit<int(vClusterDigitCollection.size()); jdigit++ ){
           RecoClusterDigit* cdigit = (RecoClusterDigit*)(vClusterDigitCollection.at(jdigit));
           RecoDigit* recodigit = (RecoDigit*)(cdigit->GetRecoDigit());
-          cluster->AddDigit(recodigit);        
+          cluster.AddDigit(recodigit);        
         }
       }
     }
@@ -738,3 +740,32 @@ std::vector<RecoDigit*>* HitCleaner::FilterByTruthInfo(std::vector<RecoDigit*>*
   
   return fFilterByTruthInfo;
 }
+
+void HitCleaner::CBCheck(std::vector<RecoDigit*>* unfilteredDigits, std::vector<RecoDigit*>* filteredDigits) {
+    //calculate unfiltered CB
+    double total_Q = 0;
+    double total_QSquared = 0;
+    for (int i=0;i<unfilteredDigits->size();i++) {
+        //if(unfilteredDigits->at(i)->GetDigitType()==RecoDigit::PMT8inch){
+        double tube_charge = unfilteredDigits->at(i)->GetCalCharge();
+        total_Q += tube_charge;
+        total_QSquared += (tube_charge * tube_charge);
+        //}
+    }
+    //FIXME: Need a method to have the 123 be equal to the number of operating detectors
+    double ucharge_balance = sqrt((total_QSquared) / (total_Q * total_Q) - (1. / 123.));
+    if (verbosity > 4) std::cout << "HitCleaner Tool: Unfiltered CB: " << ucharge_balance << std::endl;
+
+    total_Q = 0;
+    total_QSquared = 0;
+    for (int i = 0; i < filteredDigits->size(); i++) {
+        //if (filteredDigits->at(i)->GetDigitType() == RecoDigit::PMT8inch) {
+        double tube_charge = filteredDigits->at(i)->GetCalCharge();
+        total_Q += tube_charge;
+        total_QSquared += (tube_charge * tube_charge);
+        //}
+    }
+    double fcharge_balance = sqrt((total_QSquared) / (total_Q * total_Q) - (1. / 123.));
+    if (verbosity > 4) std::cout << "HitCleaner Tool: filtered CB: " << fcharge_balance << std::endl;
+
+}
diff --git a/UserTools/HitCleaner/HitCleaner.h b/UserTools/HitCleaner/HitCleaner.h
index e20cef7a0..019459d76 100644
--- a/UserTools/HitCleaner/HitCleaner.h
+++ b/UserTools/HitCleaner/HitCleaner.h
@@ -67,11 +67,12 @@ class HitCleaner: public Tool {
   std::vector<RecoDigit*>* FilterByNeighbours(std::vector<RecoDigit*>* digitlist);
   std::vector<RecoDigit*>* FilterByClusters(std::vector<RecoDigit*>* digitlist);
   std::vector<RecoDigit*>* FilterByTruthInfo(std::vector<RecoDigit*>* digitlist); //use truth information. Only for testing the code
-  std::vector<RecoCluster*>* RecoClusters(std::vector<RecoDigit*>* digitlist);
+  std::vector<RecoCluster>* RecoClusters(std::vector<RecoDigit*>* digitlist);
 
 
  private:
   void Reset();
+  void CBCheck(std::vector<RecoDigit*>* unfilteredDigits, std::vector<RecoDigit*>* filteredDigits);
   
   // running mode
   int fConfig;
@@ -120,10 +121,10 @@ class HitCleaner: public Tool {
   std::vector<RecoDigit*>* fFilterByTruthInfo;
 
   // vectors of clusters
-  std::vector<RecoCluster*>* fClusterList;
+  std::vector<RecoCluster>* fClusterList;
  
   // vector of clusters (accessible to the CStore)
-  std::vector<RecoCluster*>* fHitCleaningClusters = nullptr;  
+  std::vector<RecoCluster>* fHitCleaningClusters = nullptr;  
  	
   // true vertex
   RecoVertex* fTrueVertex = 0; 
diff --git a/UserTools/NeutronCheck/NeutronCheck.cpp b/UserTools/NeutronCheck/NeutronCheck.cpp
new file mode 100755
index 000000000..443ad40cf
--- /dev/null
+++ b/UserTools/NeutronCheck/NeutronCheck.cpp
@@ -0,0 +1,489 @@
+#include "NeutronCheck.h"
+
+NeutronCheck::NeutronCheck():Tool(){}
+
+
+bool NeutronCheck::Initialise(std::string configfile, DataModel &data){
+
+  /////////////////// Useful header ///////////////////////
+  if(configfile!="") m_variables.Initialise(configfile); // loading config file
+  //m_variables.Print();
+
+  m_variables.Get("verbosity",verbosity);
+  m_variables.Get("outfile",outfile);
+  m_variables.Get("UseClean",useClean);
+  m_variables.Get("FinderCompare",fFinderCompare);
+  m_variables.Get("ParticleInfo",fParticleInfo);
+
+  outfile+=".root";
+
+  fOutput_tfile = new TFile(outfile.c_str(), "recreate");
+  NeutCheckTree = new TTree("NeutCheckTree", "Neutron Check Tree");
+
+  NeutCheckTree->Branch("eventNumber", &fMCEventNum, "eventNumber/I");
+  
+  NeutCheckTree->Branch("eventCutStatus", &fEventStatusFlagged,"eventStatusFlagged/I");
+  NeutCheckTree->Branch("ClusterNumber",&fClusterNum);
+  NeutCheckTree->Branch("TrueEnu",&true_Enu,"true_Enu/D");
+  NeutCheckTree->Branch("TrueEmu",&true_Emu,"true_Emu/D");
+  NeutCheckTree->Branch("TrueQ2",&TrueQ2,"true_Q2/D");
+  NeutCheckTree->Branch("ClusterCount", &fClusterCount,"clusterCount/I");
+  NeutCheckTree->Branch("RecoClusters",&fRecoClusters,"RecoClusters");
+
+  NeutCheckTree->Branch("TrueNeutronMult",&fTrueNeutronMult,"trueNeutronMult/I");
+  NeutCheckTree->Branch("TrueNeutronDelayed",&fTrueNeutronDelayed,"trueNeutronDel/I");
+  NeutCheckTree->Branch("TrueNeutCapT",&fMCNeutCapTimes);
+ 
+
+  NeutCheckTree->Branch("ClusterNeutronMult",&fNeutronMult,"neutronMult/I");
+
+  NeutCheckTree->Branch("ClusterMode",&fClusterMode);
+  NeutCheckTree->Branch("ClusterPDG",&fClusterPDG);
+  NeutCheckTree->Branch("ClusterParentPDG",&fClusterParentPDG);
+  NeutCheckTree->Branch("ClusterParticleEnergy",&fClusterParticleEnergy);
+  NeutCheckTree->Branch("ClusterHits", &fClusterHits);
+  NeutCheckTree->Branch("ClusterTime",&fClusterTime);
+  NeutCheckTree->Branch("ClusterCharge",&fClusterCharge);
+  NeutCheckTree->Branch("ClusterPurity",&fClusterPurity);
+  NeutCheckTree->Branch("ClusterCB",&fClusterCB);
+  NeutCheckTree->Branch("ClusterAS0",&fClusterAS0);
+  NeutCheckTree->Branch("ClusterAS1",&fClusterAS1);
+  NeutCheckTree->Branch("ClusterAS2",&fClusterAS2);
+  //NeutCheckTree->Branch("ClusterSA",&fClusterSA);
+  NeutCheckTree->Branch("ClusterCVX",&fClusterCVX);
+  NeutCheckTree->Branch("ClusterCVY", &fClusterCVY);
+  NeutCheckTree->Branch("ClusterCVZ", &fClusterCVZ);
+  NeutCheckTree->Branch("ClusterCVR", &fClusterCVR);
+  NeutCheckTree->Branch("ClusterAMD", &fClusterAMD);
+
+
+  if(fFinderCompare){
+  NeutCheckTree->Branch("FinderClusterNumber",&fFinderClusterNum);
+  NeutCheckTree->Branch("FinderClusterCount",&fFinderClusterCount,"FinderCount/I");
+  NeutCheckTree->Branch("FinderPDG",&fFinderPDG);
+  NeutCheckTree->Branch("FinderParentPDG",&fFinderParentPDG);
+  NeutCheckTree->Branch("FinderCharge",&fFinderCharge);
+  NeutCheckTree->Branch("FinderPurity",&fFinderPurity);
+  NeutCheckTree->Branch("FinderTime",&fFinderTime);
+  }
+
+  if (fParticleInfo) {
+      //Log("Particle Information to be added.",v_debug,verbosity);
+      NeutCheckTree->Branch("ParticlePDG",&fParticlePDG);
+      NeutCheckTree->Branch("ParticleParent",&fParticleParent);
+      NeutCheckTree->Branch("ParticleStartEnergy",&fParticleStartEnergy);
+      NeutCheckTree->Branch("ParticleStartTime", &fParticleStartTime);
+      NeutCheckTree->Branch("TrueNeutCapX", &fMCNeutCapX);
+      NeutCheckTree->Branch("TrueNeutCapY", &fMCNeutCapY);
+      NeutCheckTree->Branch("TrueNeutCapZ", &fMCNeutCapZ);
+  }
+
+ 
+  
+
+  m_data= &data; //assigning transient data pointer
+  /////////////////////////////////////////////////////////////////
+
+  return true;
+}
+
+
+bool NeutronCheck::Execute(){
+
+    Log("NeutronCheck Tool Executing:",v_message,verbosity);
+
+    ResetVariables();
+
+    // MC entry number
+    m_data->Stores.at("ANNIEEvent")->Get("MCEventNum", fMCEventNum);
+
+    // MC trigger number
+    m_data->Stores.at("ANNIEEvent")->Get("MCTriggernum", fMCTriggerNum);
+
+    // ANNIE Event number
+    m_data->Stores.at("ANNIEEvent")->Get("EventNumber", fEventNumber);
+
+    auto get_flags = m_data->Stores.at("RecoEvent")->Get("EventFlagged", fEventStatusFlagged);
+
+    if ((fEventStatusFlagged) != 0 && useClean) {
+        //  if (!fEventCutStatus){
+        Log("NeutronCheck Tool: Event was flagged with one of the active cuts.", v_debug, verbosity);
+        return true;
+    }
+ 
+
+  
+
+    
+
+    //if(EventStore=="ANNIEEvent"){
+
+        GetClusterInformation();
+
+        m_data->Stores.at("ANNIEEvent")->Get("ClusterToBestParticleID", fClusterToBestParticleID);
+        m_data->Stores.at("ANNIEEvent")->Get("ClusterToBestParticlePDG", fClusterToBestParticlePDG);
+        m_data->Stores.at("ANNIEEvent")->Get("ClusterEfficiency", fClusterEfficiency);
+        m_data->Stores.at("ANNIEEvent")->Get("ClusterPurity", fClusterPurityMap);
+        m_data->Stores.at("ANNIEEvent")->Get("ClusterTotalCharge", fClusterTotalCharge);
+        //bool get_q2 = m_data->Stores["GenieInfo"]->Get("EventQ2", TrueQ2);
+
+        //if(get_q2) Log("True Q2 from GENIE: "+to_string(TrueQ2),v_debug,verbosity);
+        //else Log("Error finding Q2 info from GENIE",v_debug,verbosity);
+        //m_data->CStore.Get("ClausterMapMC", m_all_clusters);
+
+        int bestPartId;
+        for (int i=0;i<cluster_times.size();i++) {
+            fFinderClusterNum.push_back(i);
+            fFinderTime.push_back(cluster_times.at(i));
+            fFinderCharge.push_back(cluster_charges.at(i));
+            fFinderCB.push_back(cluster_cb.at(i));
+
+           bestPartId=fClusterToBestParticleID->at(cluster_times.at(i));
+           fFinderParentPDG.push_back(fMCParticles->at(bestPartId).GetParentPdg());
+           fFinderPDG.push_back(fClusterToBestParticlePDG->at(cluster_times.at(i)));
+           fFinderPurity.push_back(fClusterPurityMap->at(cluster_times.at(i)));
+           
+        }
+        fFinderClusterCount=fFinderClusterNum.size()+1;
+
+        fTrueNeutronMult=0;
+        fTrueNeutronDelayed=0;
+        Log("NeutCheck Tool: Scanning "+to_string(fMCParticles->size())+" MCParticles",v_debug,verbosity);
+        for (int i = 0; i < fMCParticles->size(); i++) {
+            if(fMCParticles->at(i).GetPdgCode()==2112 && fMCParticles->at(i).GetParentPdg()==0) {
+                fTrueNeutronMult++;
+                
+                if(fMCParticles->at(i).GetStopTime()>10000) fTrueNeutronDelayed++;
+            }
+            else if(fMCParticles->at(i).GetPdgCode()==2112 && fMCParticles->at(i).GetParentPdg()!=0){
+                Log("Found non-primary neutron at particle "+to_string(i)+" with parent PDG "+ to_string(fMCParticles->at(i).GetParentPdg()), v_debug, verbosity);
+            }
+            if(fParticleInfo){
+            fParticlePDG.push_back(fMCParticles->at(i).GetPdgCode());
+            fParticleParent.push_back(fMCParticles->at(i).GetParentPdg());
+            fParticleStartEnergy.push_back(fMCParticles->at(i).GetStartEnergy());
+            fParticleStartTime.push_back(fMCParticles->at(i).GetStartTime());
+            if (fMCParticles->at(i).GetPdgCode() == 2112) {
+                fMCNeutCapTimes.push_back(fMCParticles->at(i).GetStopTime());
+                fMCNeutCapX.push_back(fMCParticles->at(i).GetStopVertex().X());
+                fMCNeutCapY.push_back(fMCParticles->at(i).GetStopVertex().Y());
+                fMCNeutCapZ.push_back(fMCParticles->at(i).GetStopVertex().Z());
+            }
+            }
+        }
+
+        auto get_muonMCEnergy = m_data->Stores.at("RecoEvent")->Get("TrueMuonEnergy", true_Emu);
+        bool get_neutrino_energy = m_data->Stores["GenieInfo"]->Get("NeutrinoEnergy", true_Enu);
+
+        RecoVertex* truevtx = 0;
+        auto get_muonMC = m_data->Stores.at("RecoEvent")->Get("TrueVertex", truevtx);
+
+        //true_Emu*=1000;  //GeV->MeV to match other energies(unneeded, possibly)
+
+        double theta = truevtx->GetDirection().GetTheta();
+        double p = sqrt(pow(true_Emu,2)-pow(105.7,2));
+        TrueQ2=2*true_Enu*(true_Emu-p*cos(theta))-pow(105.7,2);
+        Log("NeutCheck TrueQ2: "+to_string(TrueQ2), v_debug, verbosity);
+
+    CSCheck();
+
+    //MapRecoCheck->Fill(fClusterToBestParticleID->size(), fRecoClusters->size());
+    //}
+
+   NeutCheckTree->Fill();
+  return true;
+}
+
+
+bool NeutronCheck::Finalise(){
+
+    Log("NeutronCheck Tool: Finalizing",v_message,verbosity);
+    fOutput_tfile->cd();
+    
+
+    NeutCheckTree->Write();
+
+    
+
+    fOutput_tfile->Close();
+
+    
+  return true;
+}
+
+
+bool NeutronCheck::GetClusterInformation() {
+
+    bool return_value = true;
+
+    bool get_cluster_idxN = m_data->Stores["RecoEvent"]->Get("ClusterIndicesNeutron", cluster_neutron);
+    if (!get_cluster_idxN) Log("NeutronCheck tool: No ClusterIndicesNeutron In RecoEvent!", v_error, verbosity);
+    bool get_cluster_tN = m_data->Stores["RecoEvent"]->Get("ClusterTimesNeutron", cluster_times_neutron);
+    if (!get_cluster_tN) Log("NeutronCheck tool: No ClusterTimesNeutron In RecoEvent!", v_error, verbosity);
+    bool get_cluster_qN = m_data->Stores["RecoEvent"]->Get("ClusterChargesNeutron", cluster_charges_neutron);
+    if (!get_cluster_qN) Log("NeutronCheck tool: No ClusterChargesNeutron In RecoEvent!", v_error, verbosity);
+    bool get_cluster_cbN = m_data->Stores["RecoEvent"]->Get("ClusterCBNeutron", cluster_cb_neutron);
+    if (!get_cluster_cbN) Log("NeutronCheck tool: No ClusterCBNeutron In RecoEvent!", v_error, verbosity);
+    bool get_cluster_t = m_data->Stores["RecoEvent"]->Get("ClusterTimes", cluster_times);
+    if (!get_cluster_t) Log("NeutronCheck tool: No ClusterTimes In RecoEvent!", v_error, verbosity);
+    bool get_cluster_q = m_data->Stores["RecoEvent"]->Get("ClusterCharges", cluster_charges);
+    if (!get_cluster_q) Log("NeutronCheck tool: No ClusterCharges In RecoEvent!", v_error, verbosity);
+    bool get_cluster_cb = m_data->Stores["RecoEvent"]->Get("ClusterCB", cluster_cb);
+    if (!get_cluster_cb) Log("NeutronCheck tool: No ClusterCB In RecoEvent!", v_error, verbosity);
+
+    bool goodMCParticles = m_data->Stores.at("ANNIEEvent")->Get("MCParticles", fMCParticles);
+    if (!goodMCParticles) {
+        std::cerr << "BackTracker: no MCParticles in the ANNIEEvent!" << endl;
+        return false;
+    }
+
+    return_value = (get_cluster_idxN && get_cluster_tN && get_cluster_qN && get_cluster_cbN && get_cluster_t && get_cluster_q && get_cluster_cb);
+
+    
+   
+    reco_Clusters = cluster_cb.size();
+
+
+
+    return return_value;
+
+}
+
+double NeutronCheck::ASCheck(vector<RecoDigit>*& digits, int mode) {
+    //mode config.  0 is maximal distance; 1 is max distance from highest charge
+    
+    double max_angle = 0;
+    double max_angle2 = 0;
+    if (mode == 0||mode==2) {
+        
+        double angle;
+        Position i_position,j_position;
+        for (int i = 0; i < digits->size(); i++) {
+            if (useClean && !digits->at(i).GetFilterStatus())continue;
+            i_position = digits->at(i).GetPosition();
+            for (int j = 0; j < digits->size(); j++) {
+                if (j == i)continue;
+                if (useClean && !digits->at(j).GetFilterStatus())continue;
+                j_position = digits->at(j).GetPosition();
+
+                angle = j_position.Angle(i_position);
+                if (angle > max_angle)max_angle = angle;
+
+            }
+        }
+        cout<<"Neutcheck max_angle: "<<max_angle<<endl;
+        if(mode==0)
+        return max_angle;
+    }
+
+    if (mode == 1||mode==2) {
+        double max_charge = 0;
+        double angle;
+        
+        int max_index = 0;
+        Position max_position;
+        Position i_position;
+        for (int i = 0; i < digits->size(); i++) {
+            if (digits->at(i).GetCalCharge() > max_charge) {
+                if(useClean && !digits->at(i).GetFilterStatus())continue;
+                max_charge = digits->at(i).GetCalCharge();
+                max_index = i;
+                max_position = digits->at(i).GetPosition();
+            }
+        }
+
+
+        for (int i = 0; i < digits->size(); i++) {
+            if (i == max_index)continue;
+            if (useClean && !digits->at(i).GetFilterStatus())continue;
+            i_position = digits->at(i).GetPosition();
+            angle = i_position.Angle(max_position);
+            if (angle > max_angle2)max_angle2 = angle;
+        }
+        cout<<"Neutcheck max_angle2: "<<max_angle2<<endl;
+        if(mode==1)
+        return max_angle2;
+    }
+    double ratio = 0;
+    if(max_angle>0) ratio = max_angle2/max_angle;
+    if(mode==2)return ratio;
+    
+    return 0;
+}
+
+void NeutronCheck::CSCheck() {
+
+    std::map<int,int>* fMCParticleIndexMap;
+    
+    Log("Cluster Check!",v_debug,verbosity);
+    bool goodMCParticleIndexMap = m_data->Stores.at("ANNIEEvent")->Get("TrackId_to_MCParticleIndex", fMCParticleIndexMap);
+    if (!goodMCParticleIndexMap) {
+        std::cerr << "NeutCheck: no TrackId_to_MCParticleIndex in the ANNIEEvent!" << endl;
+        return;
+    }
+    
+
+    
+    bool cluster_status = m_data->Stores.at("RecoEvent")->Get("RecoClusters", fRecoClusters);
+    if (!cluster_status) {
+        Log("Neutcheck tool found no recoclusters.",v_debug,verbosity);
+        return;
+    }
+ 
+    int cluster_size = fRecoClusters->size();
+    if (cluster_size == 0) {
+        Log("Neutcheck tool found empty recoclusterlist.",v_debug,verbosity);
+        return;
+    }
+    Log("Neutcheck tool found this many clusters: "+to_string(cluster_size),v_debug,verbosity);
+
+    int bestParent;
+    int bestParticleID=-5;
+    int bestPDG=-5;
+    double CVX,CVY,CVZ,CVR;
+
+    for (int i = 0; i < fRecoClusters->size(); i++) {
+        fClusterNum.push_back(i);
+        fClusterCount++;
+        bestParent = fRecoClusters->at(i)->calcBestParent();
+        //bestParent = fRecoClusters->at(i)->GetBestParent();
+        Log("Check 1: Particle ID "+to_string(bestParent),v_debug,verbosity);
+        if (bestParent >= fMCParticles->size()) {
+            Log("Invalid particle ID " + to_string(bestParent) + " vs number of particles: " + to_string(fMCParticles->size()) + "; skipping.",v_error,verbosity);
+            continue;
+        }
+        
+        //for (std::pair<int, int> apair : *fMCParticleIndexMap) {
+            //Log("Testing match for particle: "+to_string(apair.first) + " corresponding to ID " + to_string(apair.second), v_debug, verbosity);
+            //if(apair.second==bestParent){
+                //bestParticleID=apair.first;
+                bestPDG=fMCParticles->at(bestParent).GetPdgCode();
+                //if (bestPDG == 22) {
+
+                if (fMCParticles->at(bestParent).GetFlag() != 0) {
+                    Log("Flagged particle!  PDG "+to_string(bestPDG)+" excluding.",v_message,verbosity);
+                    bestPDG=-5;
+                }
+                    if (fMCParticles->at(bestParent).GetParentPdg() != 0) {
+                        
+                    Log("NeutCheck: PDG "+to_string(bestPDG)+" Finding parent.",v_message,verbosity);
+                    fClusterParentPDG.push_back(fMCParticles->at(bestParent).GetParentPdg());
+                    Log("Parent PDG "+to_string(fClusterParentPDG.at(fClusterParentPDG.size()-1)),v_message,verbosity);
+                    }
+                    else {
+                        Log("NeutCheck tool: no parent found.  Treating as primary. PDG "+to_string(bestPDG), v_message, verbosity);
+                        fClusterParentPDG.push_back(0);
+                        if (bestPDG == 2212 && fRecoClusters->at(i)->GetTime()>10000) {
+
+                        Log("Primary proton found ID " + to_string(bestParent) + " at time " +to_string(fRecoClusters->at(i)->GetTime())+" but particle start, stop "+to_string(fMCParticles->at(bestParent).GetStartTime())+", "+to_string(fMCParticles->at(bestParent).GetStopTime()),v_debug,verbosity);
+                        }
+                        if (bestPDG == 13 && fRecoClusters->at(i)->GetTime() > 10000) {
+
+                            Log("Primary muon found ID " + to_string(bestParent) + " at time " + to_string(fRecoClusters->at(i)->GetTime()) + " but particle start, stop " + to_string(fMCParticles->at(bestParent).GetStartTime()) + ", " + to_string(fMCParticles->at(bestParent).GetStopTime()), v_debug, verbosity);
+                        }
+                        if (bestPDG == 2112 && fRecoClusters->at(i)->GetTime() > 10000) {
+
+                            Log("Primary neutron found ID " + to_string(bestParent) + " at time " + to_string(fRecoClusters->at(i)->GetTime()) + " but particle start, stop " + to_string(fMCParticles->at(bestParent).GetStartTime()) + ", " + to_string(fMCParticles->at(bestParent).GetStopTime()), v_debug, verbosity);
+                        }
+                    }
+                    
+                //}
+                fRecoClusters->at(i)->SetPDG(bestPDG);
+                fClusterParticleEnergy.push_back(fMCParticles->at(bestParent).GetStartEnergy());
+        
+            //}
+            
+        //}
+        /*if (matchCheck == false) {
+            Log("NeutCheck: Failed to match Particle. IndexMap size: "+to_string(fMCParticleIndexMap->size()), v_error, verbosity);
+        //continue;
+        }*/
+        
+
+        //ClusterSearchCheck->Fill(fRecoClusters->at(i)->GetClusterMode());
+        //recoPDGHist->Fill(fRecoClusters->at(i)->GetPDG());
+        fClusterMode.push_back(fRecoClusters->at(i)->GetClusterMode());
+        fClusterPDG.push_back(fRecoClusters->at(i)->GetPDG());
+        if (fClusterPDG.at(fClusterPDG.size() - 1) == 2112) {
+            fNeutronMult++;
+
+        }
+
+        fClusterCharge.push_back(fRecoClusters->at(i)->GetCharge());
+        fClusterPurity.push_back(fRecoClusters->at(i)->Purity());
+        fClusterTime.push_back(fRecoClusters->at(i)->GetTime());
+        fClusterCB.push_back(fRecoClusters->at(i)->GetCB());
+        fClusterHits.push_back(fRecoClusters->at(i)->GetNDigits());
+        //if(fRecoClusters->at(i)->GetTime()>10000) recoPDGHistDelayed->Fill(fRecoClusters->at(i)->GetPDG());
+
+        //AS0RecoCheck->Fill(fRecoClusters->at(i)->GetAS(0));
+        //AS1RecoCheck->Fill(fRecoClusters->at(i)->GetAS(1));
+
+        fClusterAS0.push_back(fRecoClusters->at(i)->GetAS(0));
+        fClusterAS1.push_back(fRecoClusters->at(i)->GetAS(1));
+        fClusterAMD.push_back(fRecoClusters->at(i)->GetAMD());
+
+        if(fRecoClusters->at(i)->GetAS(0)>0){
+            //AS2RecoCheck->Fill(fRecoClusters->at(i)->GetAS(1)/fRecoClusters->at(i)->GetAS(0));
+            fClusterAS2.push_back(fClusterAS1.at(fClusterAS1.size()-1)/fClusterAS0.at(fClusterAS0.size()-1));
+        }
+        fClusterASC.push_back(fRecoClusters->at(i)->GetASC());
+        CVX= fRecoClusters->at(i)->GetCV().X();
+        CVY = fRecoClusters->at(i)->GetCV().Y();
+        CVZ = fRecoClusters->at(i)->GetCV().Z();
+        CVR = pow(pow(CVX,2)+pow(CVY,2)+pow(CVZ,2),0.5);
+        fClusterCVX.push_back(CVX);
+        fClusterCVY.push_back(CVY);
+        fClusterCVZ.push_back(CVZ);
+        fClusterCVR.push_back(CVR);
+        //fClusterSA.push_back(fRecoClusters->at(i)->GetSA());
+        
+    
+    }
+
+    return;
+}
+
+void NeutronCheck::ResetVariables() {
+    fMCEventNum=-9999;
+    fClusterNum.clear();
+    fClusterCount=0;
+    fClusterMode.clear();
+    fClusterPDG.clear();
+    fClusterParentPDG.clear();
+    fClusterParticleEnergy.clear();
+    fClusterCharge.clear();
+    fClusterPurity.clear();
+    fClusterCB.clear();
+    fClusterTime.clear();
+    fClusterAS0.clear();
+    fClusterAS1.clear();
+    fClusterAS2.clear();
+    fClusterASC.clear();
+    //fClusterSA.clear();
+    fClusterCVX.clear();
+    fClusterCVY.clear();
+    fClusterCVZ.clear();
+    fClusterCVR.clear();
+    fClusterAMD.clear();
+    fNeutronMult=0;
+    fMCNeutCapTimes.clear();
+    fMCNeutCapX.clear();
+    fMCNeutCapY.clear();
+    fMCNeutCapZ.clear();
+
+    fParticlePDG.clear();
+    fParticleParent.clear();
+    fParticleStartEnergy.clear();
+    fParticleStartTime.clear();
+
+
+    fFinderClusterNum.clear();
+    fFinderClusterCount=0;
+    fFinderPDG.clear();
+    fFinderParentPDG.clear();
+    fFinderCharge.clear();
+    fFinderPurity.clear();
+    fFinderCB.clear();
+    fFinderTime.clear();
+}
\ No newline at end of file
diff --git a/UserTools/NeutronCheck/NeutronCheck.h b/UserTools/NeutronCheck/NeutronCheck.h
new file mode 100755
index 000000000..16a01af72
--- /dev/null
+++ b/UserTools/NeutronCheck/NeutronCheck.h
@@ -0,0 +1,176 @@
+#ifndef NeutronCheck_H
+#define NeutronCheck_H
+
+#include <string>
+#include <iostream>
+#include <cmath>
+#include <vector>
+
+#include "Tool.h"
+#include "Particle.h"
+#include "Position.h"
+#include "RecoVertex.h"
+#include "RecoCluster.h"
+
+#include "TTree.h"
+#include "TFile.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TH3.h"
+#include "TGraph.h"
+#include "TGraphErrors.h"
+#include "TMath.h"
+
+
+
+/**
+ * \class NeutronCheck
+ *
+ * This is a blank template for a Tool used by the script to generate a new custom tool. Please fill out the description and author information.
+*
+* $Author: F. A. Lemmons $
+* $Date: 2025/02/28 10:44:00 $
+* Contact: franklin.lemmons@mines.sdsmt.edu
+*/
+class NeutronCheck: public Tool {
+
+
+ public:
+
+  NeutronCheck(); ///< Simple constructor
+  bool Initialise(std::string configfile,DataModel &data); ///< Initialise Function for setting up Tool resources. @param configfile The path and name of the dynamic configuration file to read in. @param data A reference to the transient data class used to pass information between Tools.
+  bool Execute(); ///< Execute function used to perform Tool purpose.
+  bool Finalise(); ///< Finalise function used to clean up resources.
+
+
+ private:
+
+	 bool GetClusterInformation();
+	 double ASCheck(std::vector<RecoDigit>*& digits,int mode);
+	 void CSCheck();
+	 void ResetVariables();
+
+	 int verbosity;
+	 std::string outfile;
+	 bool useClean=0;
+	 std::string EventStore;
+	 bool fFinderCompare;
+	 bool fParticleInfo;
+
+	 TTree* NeutCheckTree = nullptr;
+
+	 /*
+	 -EventNumber
+	 -Clusternumber
+	 -number of clusters
+	 -neutron multiplicity
+
+
+	 -Cluster PDG
+	 -Cluster Charge
+	 -Cluster Time
+	 -CB
+	 -AS 0,1,2
+
+	 */
+
+	 std::vector<int> cluster_neutron;
+	 std::vector<double> cluster_times_neutron;
+	 std::vector<double> cluster_charges_neutron;
+	 std::vector<double> cluster_cb_neutron;
+	 std::vector<double> cluster_times;
+	 std::vector<double> cluster_charges;
+	 std::vector<double> cluster_cb;
+
+	 std::map<double, int>* fClusterToBestParticleID = nullptr;
+	 std::map<double, int>* fClusterToBestParticlePDG = nullptr;
+	 std::map<double, double>* fClusterEfficiency = nullptr;
+	 std::map<double, double>* fClusterPurityMap = nullptr;
+	 std::map<double, double>* fClusterTotalCharge = nullptr;
+	 std::vector<RecoDigit>* fDigitList=nullptr;
+	 vector<RecoCluster*>* fRecoClusters;
+
+	 std::vector<MCParticle>* fMCParticles;
+	 std::vector<double> fMCNeutCapTimes;
+	 std::vector<double> fMCNeutCapX;
+	 std::vector<double> fMCNeutCapY;
+	 std::vector<double> fMCNeutCapZ;
+
+
+	 std::map<double,std::vector<MCHit>>* m_all_clusters;
+
+	 Geometry* fGeometry = nullptr;    ///< ANNIE Geometry
+
+	 /// \brief MC entry number
+	 uint32_t fMCEventNum;
+	 vector<int> fClusterNum;
+	 int fClusterCount;
+
+	 vector<int> fClusterPDG;
+	 vector<int> fClusterParentPDG;
+	 vector<double> fClusterParticleEnergy;
+	 vector<double> fClusterCharge;
+	 vector<double> fClusterPurity;
+	 vector<double> fClusterCB;
+	 vector<double> fClusterTime;
+	 vector<double> fClusterAS0;
+	 vector<double> fClusterAS1;
+	 vector<double> fClusterAS2;
+	 vector<double> fClusterASC;
+	 vector<double> fClusterAMD;
+	 //vector<double> fClusterSA;
+	 //vector<Position> fClusterCV;
+	 vector<double> fClusterCVX,fClusterCVY,fClusterCVZ,fClusterCVR;
+
+	 vector<int> fClusterMode;
+	 vector<int> fClusterHits;
+
+	 vector<int> fParticlePDG;
+	 vector<int> fParticleParent;
+	 vector<double> fParticleStartEnergy;
+	 vector<double> fParticleStartTime;
+
+	 vector<int> fFinderClusterNum;
+	 int fFinderClusterCount;
+	 vector<int> fFinderPDG;
+	 vector<int> fFinderParentPDG;
+	 vector<double> fFinderCharge;
+	 vector<double> fFinderPurity;
+	 vector<double> fFinderCB;
+	 vector<double> fFinderTime;
+
+	 
+	 int fTrueNeutronMult;
+	 int fTrueNeutronDelayed;
+	 int fNeutronMult;
+	 double true_Emu;
+	 double true_Enu;
+	 double TrueQ2;
+
+	 /// \brief trigger number
+	 uint16_t fMCTriggerNum;
+
+	 /// \brief ANNIE event number
+	 uint32_t fEventNumber;
+
+	 int fEventStatusApplied;
+	 int fEventStatusFlagged;
+
+	 
+
+	 int reco_Clusters;
+
+	 TFile* fOutput_tfile;
+
+	 //verbosity variables
+	 int v_error = 0;
+	 int v_warning = 1;
+	 int v_message = 2;
+	 int v_debug = 3;
+	 int vv_debug = 4;
+
+
+};
+
+
+#endif
diff --git a/UserTools/NeutronCheck/README.md b/UserTools/NeutronCheck/README.md
new file mode 100644
index 000000000..d2be28214
--- /dev/null
+++ b/UserTools/NeutronCheck/README.md
@@ -0,0 +1,71 @@
+# NeutronCheck
+
+NeutronCheck
+Output tool for RecoCluster information and parameters.  Creates NeutCheckTree in the outfile.root, for use making analysis and comparison plots.
+
+## Data
+
+
+From ANNIEEvent store
+**MCEventNum**
+*simple event number
+
+**ClusterToBestParticleID**
+**ClusterToBestParticlePDG**
+**ClusterEfficiency**
+**ClusterPurity**
+**ClusterTotalCharge**
+**ClusterMapMC**
+*ClusterFinder-made cluster maps from BackTracker tool for comparison
+
+**MCParticles**
+*vector of all MC Particles in the event
+
+**TrueMuonEnergy**
+*True energy of the muon for inclusion in the output
+
+**TrackId_to_MCParticleIndex**	DEPRECATED
+*map of track ID to MCParticle Index.  No longer used, as MCParticle Index is used directly in all relevant cases.
+
+
+
+From RecoEvent store
+**EventFlagged**
+*Which event selection cuts were triggered by the current event.
+
+**ClusterIndiccesNeutron**
+**ClusterTimesNeutron**
+**ClusterChargesNeutron**
+**ClusterCBNeutron**
+*ClusterFinder-made cluster maps of neutrons specifically for comparison
+
+**ClusterTimes**
+**ClusterCharges**
+**ClusterCB**
+*ClusterFinder-made cluster maps of all particles for comparison
+
+**RecoClusters**
+*vector list of all RecoClusters generated from all iterations of ClusterSearcher
+
+
+
+From GenieInfo Store
+**NeutrinoEnergy
+*True neutrino energy for inclusion in the output
+
+Output		.root file
+**NeutCheckTree**
+*Tree written into output file
+
+
+## Configuration
+
+Describe any configuration variables for NeutronCheck.
+
+```
+verbosity 1
+outfile /path/to/output/file	#'.root' is added to the end.
+UseClean 0	#whether to exclude events that do not pass event selection
+FinderCompare 0 #whether to add similar plots from ClusterFinder map-clusters to output for comparison purposes (not all parameters are written into Finder clusters)
+ParticleInfo 1	#whether to add true particle information directly from MCParticles list to output
+```
diff --git a/UserTools/Unity.h b/UserTools/Unity.h
index ece5408f4..602e4efa8 100644
--- a/UserTools/Unity.h
+++ b/UserTools/Unity.h
@@ -187,4 +187,6 @@
 #include "LAPPDLoadStore.h"
 #include "PMTWaveformSim.h"
 #include "LAPPDWaveformDisplay.h"
+#include "ClusterSearcher.h"
+#include "NeutronCheck.h"
 #include "PrintADCTraces.h"
diff --git a/configfiles/ClusterSearcher/ClusterSearcher2Config b/configfiles/ClusterSearcher/ClusterSearcher2Config
new file mode 100755
index 000000000..20af28b1b
--- /dev/null
+++ b/configfiles/ClusterSearcher/ClusterSearcher2Config
@@ -0,0 +1,22 @@
+# ClusterSearcher2 config file
+
+verbosity                  5 
+ClusterMode		     2			#searching track-like clusters	
+Config                      3			#config type: 1 = pulse height cut, 2 = +neighbour cut, 3 = +cluster cut							
+IsMC 1
+PmtMinPulseHeight           20		#minimum pulse height
+PmtNeighbourRadius          600		#digit neighbouring distance [cm]
+PmtMinNeighbourDigits       2			#minimum neighbour digits
+PmtClusterRadius            600		#digit clustering distance [cm]      
+PmtTimeWindowN              100		#neighbouring time window [ns]    
+PmtTimeWindowC              100		#clustering time window [ns]           
+PmtMinHitsPerCluster        1			#number of hits per cluster								                         
+LappdMinPulseHeight         0			#minimum pulse height                             
+LappdNeighbourRadius        250              #digit neighbouring distance [cm]              
+LappdMinNeighbourDigits     20              #minimum neighbour digits                      
+LappdClusterRadius          250              #digit clustering distance [cm]                                    
+LappdTimeWindowN            1               #neighbouring time window [ns]                 
+LappdTimeWindowC            1               #clustering time window [ns]                   
+LappdMinHitsPerCluster      5			#number of digits per cluster	
+#MinClusterDigits           50		#minimum clustered digits	(LAPPD+PMT)		
+MinClusterDigits            1			#minimum clustered digits	(PMT-only)	
diff --git a/configfiles/ClusterSearcher/ClusterSearcherConfig b/configfiles/ClusterSearcher/ClusterSearcherConfig
new file mode 100755
index 000000000..22020dbbf
--- /dev/null
+++ b/configfiles/ClusterSearcher/ClusterSearcherConfig
@@ -0,0 +1,22 @@
+# ClusterSearcher config file
+
+verbosity                   5
+ClusterMode		     1			#searching track-like clusters	
+Config                      3			#config type: 1 = pulse height cut, 2 = +neighbour cut, 3 = +cluster cut							
+IsMC 1
+PmtMinPulseHeight           20		#minimum pulse height
+PmtNeighbourRadius          60		#digit neighbouring distance [cm]
+PmtMinNeighbourDigits       2			#minimum neighbour digits
+PmtClusterRadius            60		#digit clustering distance [cm]      
+PmtTimeWindowN              10		#neighbouring time window [ns]    
+PmtTimeWindowC              100		#clustering time window [ns]           
+PmtMinHitsPerCluster        4			#number of hits per cluster								                         
+LappdMinPulseHeight         0			#minimum pulse height                             
+LappdNeighbourRadius        25              #digit neighbouring distance [cm]              
+LappdMinNeighbourDigits     20              #minimum neighbour digits                      
+LappdClusterRadius          25              #digit clustering distance [cm]                                    
+LappdTimeWindowN            1               #neighbouring time window [ns]                 
+LappdTimeWindowC            1               #clustering time window [ns]                   
+LappdMinHitsPerCluster      5			#number of digits per cluster	
+#MinClusterDigits           50		#minimum clustered digits	(LAPPD+PMT)		
+MinClusterDigits            4			#minimum clustered digits	(PMT-only)	
diff --git a/configfiles/ClusterSearcher/DigitBuilderConfig b/configfiles/ClusterSearcher/DigitBuilderConfig
new file mode 100755
index 000000000..a2aac4e86
--- /dev/null
+++ b/configfiles/ClusterSearcher/DigitBuilderConfig
@@ -0,0 +1,12 @@
+# DigitBuilder config file
+
+verbosity 5
+ParametricModel 1
+#Reading in MC files
+IsMC 1
+DigitChargeThr 20 #need to change the name
+# There are three configurations: "PMT_only", "LAPPD_only", "All"
+PhotoDetectorConfiguration PMT_only
+#File must be in /pnfs/ space when loading on Fermilab cluster
+LAPPDIDFile ./configfiles/VertexReco/PhaseIIRecoFirinne/LAPPDIDs.txt
+StripHit 0
diff --git a/configfiles/ClusterSearcher/EventSelectorConfig b/configfiles/ClusterSearcher/EventSelectorConfig
new file mode 100755
index 000000000..6fd425abd
--- /dev/null
+++ b/configfiles/ClusterSearcher/EventSelectorConfig
@@ -0,0 +1,18 @@
+# EventSelector config file
+
+verbosity 0
+MCPMTVolCut 0
+MCFVCut 1
+MCMRDCut 1
+MCPiKCut 0
+MCIsMuonCut 1
+MRDRecoCut 0
+RecoPMTVolCut 0
+RecoFVCut 0 
+ArgonFV 0
+MCIsSingleRingCut 0
+NHitCut 1
+LAPPDMult 0
+TriggerWord 0
+PromptTrigOnly 0
+SaveStatusToStore 1
diff --git a/configfiles/ClusterSearcher/LoadGenieEventConfig b/configfiles/ClusterSearcher/LoadGenieEventConfig
new file mode 100644
index 000000000..26c9051d8
--- /dev/null
+++ b/configfiles/ClusterSearcher/LoadGenieEventConfig
@@ -0,0 +1,12 @@
+verbosity 0
+FluxVersion 0  # use 0 to load genie files based on bnb_annie_0000.root etc files
+               # use 1 to load files based on beammc_annie_0000.root etc files
+#FileDir NA     # specify "NA" for newer files: full path is saved in WCSim
+#FileDir /pnfs/annie/persistent/users/vfischer/genie_files/BNB_Water_10k_22-05-17
+#FileDir /pnfs/annie/persistent/users/moflaher/genie/BNB_World_10k_11-03-18_gsimpleflux
+FileDir /Data/Neutronsim
+FilePattern gntp.1078.ghep.root
+#FilePattern LoadWCSimTool      ## use this pattern to load corresponding genie info with the LoadWCSimTool
+                               ## N.B: FileDir must still be specified for now! (WCSim files do not record their directory)
+ManualFileMatching 0
+EventOffset 0			#9900
diff --git a/configfiles/ClusterSearcher/LoadGeometryConfig b/configfiles/ClusterSearcher/LoadGeometryConfig
new file mode 100644
index 000000000..3d8ce96bc
--- /dev/null
+++ b/configfiles/ClusterSearcher/LoadGeometryConfig
@@ -0,0 +1,9 @@
+verbosity 0
+LAPPDChannelCount 60
+FACCMRDGeoFile ./configfiles/LoadGeometry/FullMRDGeometry_09_29_20.csv 
+#FACCMRDGeoFile ./configfiles/LoadGeometry/FullMRDGeometry.csv
+DetectorGeoFile ./configfiles/LoadGeometry/DetectorGeometrySpecs.csv
+LAPPDGeoFile ./configfiles/LoadGeometry/LAPPDGeometry.csv
+TankPMTGeoFile ./configfiles/LoadGeometry/FullTankPMTGeometry.csv
+TankPMTGainFile ./configfiles/LoadGeometry/ChannelSPEGains_BeamRun20192020.csv
+AuxiliaryChannelFile ./configfiles/LoadGeometry/AuxChannels.csv
diff --git a/configfiles/ClusterSearcher/LoadWCSimConfig b/configfiles/ClusterSearcher/LoadWCSimConfig
new file mode 100644
index 000000000..9bd2001a0
--- /dev/null
+++ b/configfiles/ClusterSearcher/LoadWCSimConfig
@@ -0,0 +1,17 @@
+#LoadWCSim Config File
+# all variables retrieved with m_variables.Get() must be defined here!
+verbose 0
+
+#InputFile /Data/Ambesim/pmt/*.root
+#InputFile /Data/HighEnergy/wcsim-thousand-high_0.root
+InputFile /Data/0824Beamsim/pmt/*.root
+
+WCSimVersion 3               ## should reflect the WCSim version of the files being loaded
+HistoricTriggeroffset 0      ## time offset of digits relative to the trigger
+UseDigitSmearedTime 1        ## whether to use smeared digit time (T), or true time of first photon (F)
+LappdNumStrips 60            ## num channels to construct from each LAPPD
+LappdStripLength 100         ## relative x position of each LAPPD strip, for dual-sided readout [mm]
+LappdStripSeparation 10      ## stripline separation, for calculating relative y position of each LAPPD strip [mm]
+PMTMask configfiles/BeamClusterAnalysisMC/DeadPMTIDs_p2v7.txt ## Which PMTs should be masked out? / are dead?
+ChankeyToPMTIDMap ./configfiles/BeamClusterAnalysisMC/Chankey_WCSimID_v7.txt
+SplitSubTriggers 0		# should subtriggers be loaded in separate Execute steps?
diff --git a/configfiles/ClusterSearcher/LoadWCSimLAPPDConfig b/configfiles/ClusterSearcher/LoadWCSimLAPPDConfig
new file mode 100644
index 000000000..f39db46d8
--- /dev/null
+++ b/configfiles/ClusterSearcher/LoadWCSimLAPPDConfig
@@ -0,0 +1,11 @@
+#LoadWCSimLAPPD Config File
+# all variables retrieved with m_variables.Get() must be defined here!
+verbose 1
+
+#InputFile /Data/Ambesim/lappd/*.root
+#InputFile /Data/HighEnergy/wcsim-thousand-high_lappd_0.root
+InputFile /Data/0824Beamsim/lappd/*.root
+
+WCSimVersion 3               ## should reflect the WCSim version of the files being loaded
+InnerStructureRadius 1.3545  ## octagonal inner structure radius in m (from drawings 106.64")
+DrawDebugGraphs 0            ## whether to draw TPolyMarker3D's of hits
diff --git a/configfiles/ClusterSearcher/MCParticlePropertiesConfig b/configfiles/ClusterSearcher/MCParticlePropertiesConfig
new file mode 100644
index 000000000..714c40f32
--- /dev/null
+++ b/configfiles/ClusterSearcher/MCParticlePropertiesConfig
@@ -0,0 +1,3 @@
+# MCParticleProperties configuration file
+
+verbosity 5
diff --git a/configfiles/ClusterSearcher/MCRecoEventLoaderConfig b/configfiles/ClusterSearcher/MCRecoEventLoaderConfig
new file mode 100644
index 000000000..cba5524e7
--- /dev/null
+++ b/configfiles/ClusterSearcher/MCRecoEventLoaderConfig
@@ -0,0 +1,11 @@
+# MCRecoEventLoader config file
+
+verbosity 5
+GetPionKaonInfo 1
+GetNRings 1
+ParticleID 13
+DoParticleSelection 0
+xshift 0.0
+yshift 14.46469
+zshift -168.1
+IsMC 1
diff --git a/configfiles/ClusterSearcher/NeutronCheckConfig b/configfiles/ClusterSearcher/NeutronCheckConfig
new file mode 100644
index 000000000..0658b1a35
--- /dev/null
+++ b/configfiles/ClusterSearcher/NeutronCheckConfig
@@ -0,0 +1,6 @@
+verbosity 5
+outfile /Data/0824Beamsim/NeutronCheck_CS
+UseClean 1
+FinderCompare 0
+ParticleInfo 0
+VertexInfo 0
diff --git a/configfiles/ClusterSearcher/ToolChainConfig b/configfiles/ClusterSearcher/ToolChainConfig
new file mode 100644
index 000000000..df77bb067
--- /dev/null
+++ b/configfiles/ClusterSearcher/ToolChainConfig
@@ -0,0 +1,23 @@
+#ToolChain dynamic setup file
+
+##### Runtime Paramiters #####
+verbose 1
+error_level 0 # 0= do not exit, 1= exit on unhandeled errors only, 2= exit on unhandeled errors and handeled errors
+attempt_recover 1
+
+###### Logging #####
+log_mode Interactive # Interactive=cout , Remote= remote logging system "serservice_name Remote_Logging" , Local = local file log;
+log_local_path ./log
+log_service LogStore
+
+###### Service discovery #####
+service_publish_sec -1
+service_kick_sec -1
+
+##### Tools To Add #####
+Tools_File ./configfiles/ClusterSearcher/ToolsConfig
+
+##### Run Type #####
+Inline -1
+Interactive 0
+
diff --git a/configfiles/ClusterSearcher/ToolsConfig b/configfiles/ClusterSearcher/ToolsConfig
new file mode 100755
index 000000000..7da2cdb68
--- /dev/null
+++ b/configfiles/ClusterSearcher/ToolsConfig
@@ -0,0 +1,11 @@
+LoadGeometry LoadGeometry ./configfiles/ClusterSearcher/LoadGeometryConfig
+LoadWCSim LoadWCSim ./configfiles/ClusterSearcher/LoadWCSimConfig
+LoadWCSimLAPPD LoadWCSimLAPPD ./configfiles/ClusterSearcher/LoadWCSimLAPPDConfig
+LoadGenieEvent LoadGenieEvent ./configfiles/ClusterSearcher/LoadGenieEventConfig
+MCRecoEventLoader MCRecoEventLoader ./configfiles/ClusterSearcher/MCRecoEventLoaderConfig
+MCParticleProperties MCParticleProperties ./configfiles/ClusterSearcher/MCParticlePropertiesConfig
+DigitBuilder DigitBuilder ./configfiles/ClusterSearcher/DigitBuilderConfig
+ClusterSearcher ClusterSearcher ./configfiles/ClusterSearcher/ClusterSearcherConfig
+ClusterSearcher2 ClusterSearcher ./configfiles/ClusterSearcher/ClusterSearcher2Config
+EventSelector EventSelector ./configfiles/ClusterSearcher/EventSelectorConfig
+NeutronCheck NeutronCheck ./configfiles/ClusterSearcher/NeutronCheckConfig