diff --git a/Recon/include/Recon/ParticleFlow.h b/Recon/include/Recon/ParticleFlow.h
index dc2d2911b..b61e93bfb 100644
--- a/Recon/include/Recon/ParticleFlow.h
+++ b/Recon/include/Recon/ParticleFlow.h
@@ -13,6 +13,13 @@
 #include "Framework/EventProcessor.h"  //Needed to declare processor
 #include "TGraph.h"
 
+#include "SimCore/Event/SimParticle.h"
+#include "SimCore/Event/SimTrackerHit.h"
+#include "Recon/Event/CaloCluster.h"
+#include "Ecal/Event/EcalCluster.h"
+#include "Hcal/Event/HcalCluster.h"
+#include "Recon/Event/PFCandidate.h"
+
 namespace recon {
 
 /**
@@ -36,6 +43,10 @@ class ParticleFlow : public framework::Producer {
 
   virtual void onProcessEnd();
 
+  void FillCandTrack(ldmx::PFCandidate &cand, const ldmx::SimTrackerHit &tk);
+  void FillCandEMCalo(ldmx::PFCandidate &cand, const ldmx::CaloCluster &em);
+  void FillCandHadCalo(ldmx::PFCandidate &cand, const ldmx::CaloCluster &had);
+
  private:
   // specific verbosity of this producer
   int verbose_{0};
@@ -49,6 +60,8 @@ class ParticleFlow : public framework::Producer {
   std::string inputTrackCollName_;
   // name of collection for PF outputs
   std::string outputCollName_;
+  // configuration
+  bool singleParticle_;
 };
 }  // namespace recon
 
diff --git a/Recon/python/pfReco.py b/Recon/python/pfReco.py
index 372fd9374..30ce40f1c 100644
--- a/Recon/python/pfReco.py
+++ b/Recon/python/pfReco.py
@@ -49,6 +49,7 @@ def __init__(self, name='PFlow') :
         self.inputHcalCollName  = 'PFHcalClusters'
         self.inputTrackCollName = 'PFTracks'
         self.outputCollName     = 'PFCandidates'
+        self.singleParticle     = False
   
 class pfTruthProducer(ldmxcfg.Producer) :
     """Configuration for track selector for particle reco"""
diff --git a/Recon/src/Recon/DBScanClusterBuilder.cxx b/Recon/src/Recon/DBScanClusterBuilder.cxx
index 96ccf0004..6c33ce235 100644
--- a/Recon/src/Recon/DBScanClusterBuilder.cxx
+++ b/Recon/src/Recon/DBScanClusterBuilder.cxx
@@ -82,9 +82,9 @@ namespace recon {
   float e(0),x(0),y(0),z(0),xx(0),yy(0),zz(0),n(0);
   float w = 1; // weight
   float sumw = 0;
-  std::vector<float> xvals{};
-  std::vector<float> yvals{};
-  std::vector<float> zvals{};
+  // std::vector<float> xvals{};
+  // std::vector<float> yvals{};
+  // std::vector<float> zvals{};
   std::vector<float> raw_xvals{};
   std::vector<float> raw_yvals{};
   std::vector<float> raw_zvals{};
@@ -102,9 +102,9 @@ namespace recon {
     zz += w * h->getZPos() * h->getZPos();
     n += 1;
     sumw += w;
-    xvals.push_back(x);
-    yvals.push_back(y);
-    zvals.push_back(z);
+    // xvals.push_back(x); // unused
+    // yvals.push_back(y);
+    // zvals.push_back(z);
     raw_xvals.push_back(h->getXPos());
     raw_yvals.push_back(h->getYPos());
     raw_zvals.push_back(h->getZPos());
@@ -128,20 +128,28 @@ namespace recon {
   cl->setHitValsZ(raw_zvals);
   cl->setHitValsE(raw_evals);
 
-  if (xvals.size()>2){
-    for(int i=0; i<xvals.size();i++){ // mean subtract
-      xvals[i] = xvals[i] - x;
-      yvals[i] = yvals[i] - y;
-      zvals[i] = zvals[i] - z;
+  if (raw_xvals.size()>2){
+    // skip fits for 'vertical' clusters
+    std::vector<float> sortedZ = raw_zvals;
+    std::sort(sortedZ.begin(), sortedZ.end());
+    if (sortedZ[sortedZ.size()-1] - sortedZ[0] > 1e3){
+
+      for(int i=0; i<raw_xvals.size();i++){ // mean subtract
+	raw_xvals[i] = raw_xvals[i] - x;
+	raw_yvals[i] = raw_yvals[i] - y;
+	raw_zvals[i] = raw_zvals[i] - z;
+      }
+
+      TGraph gxz(raw_zvals.size(), raw_zvals.data(), raw_xvals.data());
+      auto r_xz = gxz.Fit("pol1","SQ"); // p0 + x*p1
+      cl->setDXDZ( r_xz->Value(1) );
+      cl->setEDXDZ( r_xz->ParError(1) );
+    
+      TGraph gyz(raw_zvals.size(), raw_zvals.data(), raw_yvals.data());
+      auto r_yz = gyz.Fit("pol1","SQ"); // p0 + x*p1
+      cl->setDYDZ( r_yz->Value(1) );
+      cl->setEDYDZ( r_yz->ParError(1) );
     }
-    TGraph gxz(zvals.size(), zvals.data(), xvals.data());
-    auto r_xz = gxz.Fit("pol1","SQ"); // p0 + x*p1
-    cl->setDXDZ( r_xz->Value(1) );
-    cl->setEDXDZ( r_xz->ParError(1) );
-    TGraph gyz(zvals.size(), zvals.data(), yvals.data());
-    auto r_yz = gyz.Fit("pol1","SQ"); // p0 + x*p1
-    cl->setDYDZ( r_yz->Value(1) );
-    cl->setEDYDZ( r_yz->ParError(1) );
   }
   return;
 }
diff --git a/Recon/src/Recon/PFTrackProducer.cxx b/Recon/src/Recon/PFTrackProducer.cxx
index ec5db2753..3d1b80965 100644
--- a/Recon/src/Recon/PFTrackProducer.cxx
+++ b/Recon/src/Recon/PFTrackProducer.cxx
@@ -10,6 +10,11 @@ void PFTrackProducer::configure(framework::config::Parameters& ps) {
   outputTrackCollName_ = ps.getParameter<std::string>("outputTrackCollName"); 
 }
 
+double getP(const ldmx::SimTrackerHit &tk){
+  std::vector<double> pxyz = tk.getMomentum();
+  return sqrt(pow(pxyz[0],2)+pow(pxyz[1],2)+pow(pxyz[2],2));
+}
+
 void PFTrackProducer::produce(framework::Event& event) {
 
   if (!event.exists(inputTrackCollName_)) return;
@@ -26,7 +31,7 @@ void PFTrackProducer::produce(framework::Event& event) {
   }
   std::sort(pfTracks.begin(), pfTracks.end(),
 	    [](ldmx::SimTrackerHit a, ldmx::SimTrackerHit b) {
-	      return a.getEnergy() > b.getEnergy();
+	      return getP(a) > getP(b);
 	    });
   event.add(outputTrackCollName_, pfTracks);
 }
diff --git a/Recon/src/Recon/ParticleFlow.cxx b/Recon/src/Recon/ParticleFlow.cxx
index 689b356c9..02c848955 100644
--- a/Recon/src/Recon/ParticleFlow.cxx
+++ b/Recon/src/Recon/ParticleFlow.cxx
@@ -1,21 +1,19 @@
 #include "Recon/ParticleFlow.h"
 
-#include "SimCore/Event/SimParticle.h"
-#include "SimCore/Event/SimTrackerHit.h"
-#include "Recon/Event/CaloCluster.h"
-#include "Ecal/Event/EcalCluster.h"
-#include "Hcal/Event/HcalCluster.h"
-#include "Recon/Event/PFCandidate.h"
 #include <vector>
 
 namespace recon {
 
 void ParticleFlow::configure(framework::config::Parameters& ps) {
+  // I/O
   inputEcalCollName_ = ps.getParameter<std::string>("inputEcalCollName");
   inputHcalCollName_ = ps.getParameter<std::string>("inputHcalCollName");
   inputTrackCollName_ = ps.getParameter<std::string>("inputTrackCollName");
   outputCollName_ = ps.getParameter<std::string>("outputCollName"); 
-  //from jason
+  // Algorithm configuration
+  singleParticle_ = ps.getParameter<bool>("singleParticle"); 
+
+  // Calibration factors, from jason
   std::vector<float> em1{250.0,750.0,1250.0,1750.0,2250.0,2750.0,3250.0,3750.0,4250.0,4750.0,5250.0,5750.};
   std::vector<float> em2{1.175,1.02,0.99,0.985,0.975,0.975,0.96,0.94,0.87,0.8,0.73,0.665};
   std::vector<float> h1{25.0,75.0,125.0,175.0,225.0,275.0,325.0,375.0,425.0};
@@ -24,6 +22,75 @@ void ParticleFlow::configure(framework::config::Parameters& ps) {
   hCorr_ = new TGraph(h1.size(), h1.data(), h2.data());
 }
 
+void ParticleFlow::FillCandTrack(ldmx::PFCandidate &cand, const ldmx::SimTrackerHit &tk){
+  // TODO: smear
+  std::vector<float> xyz = tk.getPosition();
+  std::vector<double> pxyz = tk.getMomentum();
+  float ecalZ = 248;
+  float ecalX = xyz[0] + pxyz[0]/pxyz[2] * (ecalZ - xyz[2]);
+  float ecalY = xyz[1] + pxyz[1]/pxyz[2] * (ecalZ - xyz[2]);    
+  cand.setEcalPositionXYZ(ecalX, ecalY, ecalZ);
+  cand.setTrackPxPyPz(pxyz[0], pxyz[1], pxyz[2]);
+  // also use this object to set truth info
+  cand.setTruthEcalXYZ(ecalX, ecalY, ecalZ);
+  cand.setTruthPxPyPz(pxyz[0], pxyz[1], pxyz[2]);
+  float m2 = pow(tk.getEnergy(),2) - pow(pxyz[0],2) - pow(pxyz[1],2) - pow(pxyz[2],2);
+  if(m2<0) m2=0;
+  cand.setTruthMass(sqrt(m2));
+  cand.setTruthEnergy(tk.getEnergy());
+  cand.setTruthPdgId(tk.getPdgID());
+  cand.setPID( cand.getPID() | 1 );
+}
+void ParticleFlow::FillCandEMCalo(ldmx::PFCandidate &cand, const ldmx::CaloCluster &em){
+  float corr = 1.;
+  float e = em.getEnergy();
+  if (e<eCorr_->GetX()[0]){
+    corr = eCorr_->GetX()[0];
+  } else if (e>eCorr_->GetX()[eCorr_->GetN()-1]){
+    corr = eCorr_->GetX()[eCorr_->GetN()-1];
+  } else {
+    corr = eCorr_->Eval(e);
+  }
+  cand.setEcalEnergy( e*corr );
+  cand.setEcalRawEnergy( e );
+  cand.setEcalClusterXYZ(em.getCentroidX(),
+		       em.getCentroidY(),
+		       em.getCentroidZ());
+  cand.setEcalClusterEXYZ(em.getRMSX(),
+			em.getRMSY(),
+			em.getRMSZ());
+  cand.setEcalClusterDXDZ(em.getDXDZ());
+  cand.setEcalClusterDYDZ(em.getDYDZ());
+  cand.setEcalClusterEDXDZ(em.getEDXDZ());
+  cand.setEcalClusterEDYDZ(em.getEDYDZ());
+  cand.setPID( cand.getPID() | 2 );
+}
+void ParticleFlow::FillCandHadCalo(ldmx::PFCandidate &cand, const ldmx::CaloCluster &had){
+  float corr = 1.;
+  float e = had.getEnergy();
+  if (e<hCorr_->GetX()[0]){
+    corr = hCorr_->GetX()[0];
+  } else if (e>hCorr_->GetX()[hCorr_->GetN()-1]){
+    corr = hCorr_->GetX()[hCorr_->GetN()-1];
+  } else {
+    corr = hCorr_->Eval(e);
+  }
+  cand.setHcalEnergy( e*corr );
+  cand.setHcalRawEnergy( e );
+  cand.setHcalClusterXYZ(had.getCentroidX(),
+			 had.getCentroidY(),
+			 had.getCentroidZ());
+  cand.setHcalClusterEXYZ(had.getRMSX(),
+			  had.getRMSY(),
+			  had.getRMSZ());
+  cand.setHcalClusterDXDZ(had.getDXDZ());
+  cand.setHcalClusterDYDZ(had.getDYDZ());
+  cand.setHcalClusterEDXDZ(had.getEDXDZ());
+  cand.setHcalClusterEDYDZ(had.getEDYDZ());
+  cand.setPID( cand.getPID() | 4 );
+}
+
+
 void ParticleFlow::produce(framework::Event& event) {
 
   if (!event.exists(inputTrackCollName_)) return;
@@ -32,142 +99,267 @@ void ParticleFlow::produce(framework::Event& event) {
   const auto ecalClusters = event.getCollection<ldmx::CaloCluster>(inputEcalCollName_);
   const auto hcalClusters = event.getCollection<ldmx::CaloCluster>(inputHcalCollName_);
   const auto tracks = event.getCollection<ldmx::SimTrackerHit>(inputTrackCollName_);
-
-  //std::cout << ecalClusters.size() << " " << hcalClusters.size() << " " << tracks.size() << std::endl;
-
-  /*
-    1. Build links maps at the Tk/Ecal and Hcal/Hcal interfaces
-    2. Categorize tracks as: Ecal-matched, (Side) Hcal-matched, unmatched
-    3. Categorize Ecal clusters as: Hcal-matched, unmatched
-    4a. (Upstream?) Categorize tracks with dedx?
-    4b. (Upstream?) Categorize ecal clusters as: EM/Had-like
-    4c. (Upstream?) Categorize hcal clusters as: EM/Had-like
-    5. Build candidates by category, moving from Tk-Ecal-Hcal
-   */
-
-  // track-calo linking
-  std::map<int, std::vector<int> > tkCaloMap;
-  for(int i=0; i<tracks.size(); i++){
-    const auto& tk = tracks[i];
-    for(int j=0; j<ecalClusters.size(); j++){
-      const auto& ecal = ecalClusters[j];
-      // TODO: add the matching logic here...
-      bool isMatch = true;
-      if(isMatch){
-	if (tkCaloMap.count(i)) tkCaloMap[i].push_back(j);
-	else tkCaloMap[i] = {j};
+  
+  std::vector<ldmx::PFCandidate> pfCands;
+  if(!singleParticle_){
+    /*
+      1. Build links maps at the Tk/Ecal and Hcal/Hcal interfaces
+      2. Categorize tracks as: Ecal-matched, (Side) Hcal-matched, unmatched
+      3. Categorize Ecal clusters as: Hcal-matched, unmatched
+      4a. (Upstream?) Categorize tracks with dedx?
+      4b. (Upstream?) Categorize ecal clusters as: EM/Had-like
+      4c. (Upstream?) Categorize hcal clusters as: EM/Had-like
+      5. Build candidates by category, moving from Tk-Ecal-Hcal
+    */
+    
+    //
+    // track-calo linking
+    //
+    std::map<int, std::vector<int> > tkCaloMap;
+    std::map<int, std::vector<int> > caloTkMap;
+    std::map<std::pair<int,int>, float > tkEMDist;
+    //std::vector<int> unmatchedTracks;
+    for(int i=0; i<tracks.size(); i++){
+      const auto& tk = tracks[i];
+      const std::vector<float> xyz = tk.getPosition();
+      const std::vector<double> pxyz = tk.getMomentum();
+      const float p = sqrt(pow(pxyz[0],2)+pow(pxyz[1],2)+pow(pxyz[2],2));
+      //float bestMatchVal = 9e9;
+      for(int j=0; j<ecalClusters.size(); j++){
+	const auto& ecal = ecalClusters[j];
+	// Matching logic
+	const float ecalClusZ = ecal.getCentroidZ();	
+	const float tkXAtClus = xyz[0] + pxyz[0]/pxyz[2] * (ecalClusZ - xyz[2]);
+	const float tkYAtClus = xyz[1] + pxyz[1]/pxyz[2] * (ecalClusZ - xyz[2]);
+	float dist = hypot((tkXAtClus - ecal.getCentroidX()) / std::max(1.0, ecal.getRMSX()),
+			   (tkYAtClus - ecal.getCentroidY()) / std::max(1.0, ecal.getRMSY()));
+	tkEMDist[{i,j}] = dist;
+	bool isMatch = (dist<2) && (ecal.getEnergy() > 0.3*p && ecal.getEnergy() < 2*p);
+	//if (isMatch && dist < bestMatchVal) bestMatchVal = dist;
+	if(isMatch){
+	  if (tkCaloMap.count(i)) tkCaloMap[i].push_back(j);
+	  else tkCaloMap[i] = {j};
+	  if (caloTkMap.count(j)) caloTkMap[j].push_back(i);
+	  else caloTkMap[j] = {i};
+	}
       }
     }
-  }
   
-  // em-hadcalo linking
-  std::map<int, std::vector<int> > emHadCaloMap;
-  for(int i=0; i<ecalClusters.size(); i++){
-    const auto& ecal = ecalClusters[i];
-    for(int j=0; j<hcalClusters.size(); j++){
-      const auto& hcal = hcalClusters[j];
-      // TODO: add the matching logic here...
-      bool isMatch = true;
-      if(isMatch){
-	if (emHadCaloMap.count(i)) emHadCaloMap[i].push_back(j);
-	else emHadCaloMap[i] = {j};
+    // em-hadcalo linking
+    std::map<int, std::vector<int> > emHadCaloMap;
+    std::map<std::pair<int,int>, float > EMHadDist;
+    for(int i=0; i<ecalClusters.size(); i++){
+      const auto& ecal = ecalClusters[i];
+      for(int j=0; j<hcalClusters.size(); j++){
+	const auto& hcal = hcalClusters[j];
+	// TODO: matching logic
+	const float xAtHClus = ecal.getCentroidX() + ecal.getDXDZ() * (hcal.getCentroidZ() - ecal.getCentroidZ());
+	const float yAtHClus = ecal.getCentroidY() + ecal.getDYDZ() * (hcal.getCentroidZ() - ecal.getCentroidZ());
+	float dist = sqrt( pow(xAtHClus - hcal.getCentroidX(),2) / std::max(1.0, pow(hcal.getRMSX(),2)+pow(ecal.getRMSX(),2)) +
+			   pow(yAtHClus - hcal.getCentroidY(),2) / std::max(1.0, pow(hcal.getRMSY(),2)+pow(ecal.getRMSY(),2)));
+	EMHadDist[{i,j}] = dist;
+	bool isMatch = (dist<2);
+	if(isMatch){
+	  if (emHadCaloMap.count(i)) emHadCaloMap[i].push_back(j);
+	  else emHadCaloMap[i] = {j};
+	}
       }
     }
-  }
 
-  // tk-hadcalo linking (Side HCal)
-  std::map<int, std::vector<int> > tkHadCaloMap;
-  for(int i=0; i<tracks.size(); i++){
-    const auto& tk = tracks[i];
-    for(int j=0; j<hcalClusters.size(); j++){
-      const auto& hcal = hcalClusters[j];
-      // TODO: add the matching logic here...
-      bool isMatch = true;
-      if(isMatch){
-	if (tkHadCaloMap.count(i)) tkHadCaloMap[i].push_back(j);
-	else tkHadCaloMap[i] = {j};
+    // NOT YET IMPLEMENTED...
+    // tk-hadcalo linking (Side HCal)
+    std::map<int, std::vector<int> > tkHadCaloMap;
+    // for(int i=0; i<tracks.size(); i++){
+    //   const auto& tk = tracks[i];
+    //   for(int j=0; j<hcalClusters.size(); j++){
+    // 	const auto& hcal = hcalClusters[j];
+    // 	// TODO: add the matching logic here...
+    // 	bool isMatch = true;
+    // 	if(isMatch){
+    // 	  if (tkHadCaloMap.count(i)) tkHadCaloMap[i].push_back(j);
+    // 	  else tkHadCaloMap[i] = {j};
+    // 	}
+    //   }
+    // }
+    
+    //
+    // track / ecal cluster arbitration
+    //
+    std::vector<bool> tkIsEMLinked (tracks.size(), false);
+    std::vector<bool> EMIsTkLinked (ecalClusters.size(), false);
+    std::map<int,int> tkEMPairs{};
+    for(int i=0; i<tracks.size(); i++){
+      if( tkCaloMap.count(i) ){
+	// pick first (highest-energy) unused matching cluster
+	for(int em_idx : tkCaloMap[i]){
+	  if(!EMIsTkLinked[em_idx]){
+	    EMIsTkLinked[ em_idx ] = true;
+	    tkIsEMLinked[ i ] = true;
+	    tkEMPairs[i] = em_idx;
+	    break;
+	  }
+	}
       }
     }
-  }
-  
-  std::vector<ldmx::PFCandidate> pfCands;
-  // Single-particle builder
-  ldmx::PFCandidate pf;
-  int pid=0;
-  if (tracks.size()){
-    const auto& tk = tracks[0];
-    // TODO: smear
-    std::vector<float> xyz = tk.getPosition();
-    std::vector<double> pxyz = tk.getMomentum();
-    float ecalZ = 248;
-    float ecalX = xyz[0] + pxyz[0]/pxyz[2] * (ecalZ - xyz[2]);
-    float ecalY = xyz[1] + pxyz[1]/pxyz[2] * (ecalZ - xyz[2]);    
-    pf.setEcalPositionXYZ(ecalX, ecalY, ecalZ);
-    pf.setTrackPxPyPz(pxyz[0], pxyz[1], pxyz[2]);
-    pid += 1;
-    // also use this object to set truth info
-    pf.setTruthEcalXYZ(ecalX, ecalY, ecalZ);
-    pf.setTruthPxPyPz(pxyz[0], pxyz[1], pxyz[2]);
-    float m2 = pow(tk.getEnergy(),2) - pow(pxyz[0],2) - pow(pxyz[1],2) - pow(pxyz[2],2);
-    if(m2<0) m2=0;
-    pf.setTruthMass(sqrt(m2));
-    pf.setTruthEnergy(tk.getEnergy());
-    pf.setTruthPdgId(tk.getPdgID());
-  }
-  if (ecalClusters.size()){
-    const auto& em = ecalClusters[0];
-    float corr = 1.;
-    float e = em.getEnergy();
-    if (e<eCorr_->GetX()[0]){
-      corr = eCorr_->GetX()[0];
-    } else if (e>eCorr_->GetX()[eCorr_->GetN()-1]){
-      corr = eCorr_->GetX()[eCorr_->GetN()-1];
-    } else {
-      corr = eCorr_->Eval(e);
+
+    // track / ecal cluster arbitration
+    std::vector<bool> EMIsHadLinked (ecalClusters.size(), false);
+    std::vector<bool> HadIsEMLinked (hcalClusters.size(), false);
+    std::map<int,int> EMHadPairs{};
+    for(int i=0; i<ecalClusters.size(); i++){
+      if( emHadCaloMap.count(i) ){
+	// pick first (highest-energy) unused matching cluster
+	for(int had_idx : emHadCaloMap[i]){
+	  if(!HadIsEMLinked[had_idx]){
+	    HadIsEMLinked[ had_idx ] = true;
+	    EMIsHadLinked[ i ] = true;
+	    EMHadPairs[i] = had_idx;
+	    break;
+	  }
+	}
+      }
     }
-    pf.setEcalEnergy( e*corr );
-    pf.setEcalRawEnergy( e );
-    pf.setEcalClusterXYZ(em.getCentroidX(),
-			 em.getCentroidY(),
-			 em.getCentroidZ());
-    pf.setEcalClusterEXYZ(em.getRMSX(),
-			  em.getRMSY(),
-			  em.getRMSZ());
-    pf.setEcalClusterDXDZ(em.getDXDZ());
-    pf.setEcalClusterDYDZ(em.getDYDZ());
-    pf.setEcalClusterEDXDZ(em.getEDXDZ());
-    pf.setEcalClusterEDYDZ(em.getEDYDZ());
-			 
-    pid += 2;
-  }
-  if (hcalClusters.size()){
-    const auto& had = hcalClusters[0];
-    float corr = 1.;
-    float e = had.getEnergy();
-    if (e<hCorr_->GetX()[0]){
-      corr = hCorr_->GetX()[0];
-    } else if (e>hCorr_->GetX()[hCorr_->GetN()-1]){
-      corr = hCorr_->GetX()[hCorr_->GetN()-1];
-    } else {
-      corr = hCorr_->Eval(e);
+
+    // can consider combining satellite clusters here...
+    // define some "primary cluster" ID criterion
+    //   and can add fails to the primaries
+
+    //
+    // Begin building pf candidates from tracks
+    //
+    // std::vector<ldmx::PFCandidate> chargedMatch;
+    // std::vector<ldmx::PFCandidate> chargedUnmatch;
+    for(int i=0; i<tracks.size(); i++){
+      ldmx::PFCandidate cand;
+      FillCandTrack(cand, tracks[i]);
+      if(!tkIsEMLinked[i]){
+	//chargedUnmatch.push_back(cand);
+      } else {
+	FillCandEMCalo(cand, ecalClusters[tkEMPairs[i]]);
+	if(EMIsHadLinked[tkEMPairs[i]]){
+	  FillCandHadCalo(cand, hcalClusters[EMHadPairs[tkEMPairs[i]]]);
+	}
+	//chargedMatch.push_back(cand);	
+      }
+      pfCands.push_back(cand);
     }
-    pf.setHcalEnergy( e*corr );
-    pf.setHcalRawEnergy( e );
-    pf.setHcalClusterXYZ(had.getCentroidX(),
-			 had.getCentroidY(),
-			 had.getCentroidZ());
-    pf.setHcalClusterEXYZ(had.getRMSX(),
-			  had.getRMSY(),
-			  had.getRMSZ());
-    pf.setHcalClusterDXDZ(had.getDXDZ());
-    pf.setHcalClusterDYDZ(had.getDYDZ());
-    pf.setHcalClusterEDXDZ(had.getEDXDZ());
-    pf.setHcalClusterEDYDZ(had.getEDYDZ());
-    pid += 4;
+    // std::vector<ldmx::PFCandidate> emMatch;
+    // std::vector<ldmx::PFCandidate> emUnmatch;
+    for(int i=0; i<ecalClusters.size(); i++){
+      if(EMIsTkLinked[i]) continue;
+      ldmx::PFCandidate cand;
+      FillCandEMCalo(cand, ecalClusters[i]);
+      if(EMIsHadLinked[tkEMPairs[i]]){
+	FillCandHadCalo(cand, hcalClusters[EMHadPairs[i]]);
+	//emMatch.push_back(cand);	
+      } else {
+	//emUnmatch.push_back(cand);	
+      }
+      pfCands.push_back(cand);
+
+    }
+    std::vector<ldmx::PFCandidate> hadOnly;
+    for(int i=0; i<hcalClusters.size(); i++){
+      if(HadIsEMLinked[i]) continue;
+      ldmx::PFCandidate cand;
+      FillCandHadCalo(cand, hcalClusters[i]);
+      // hadOnly.push_back(cand);
+      pfCands.push_back(cand);
+    }
+    
+
+    // // track / ecal cluster arbitration
+    // std::vector<ldmx::PFCandidate> caloMatchedTks;
+    // std::vector<ldmx::PFCandidate> unmatchedTks;
+    // std::vector<bool> emUsed (ecalClusters.size(), false);
+    // for(int i=0; i<tracks.size(); i++){
+    //   ldmx::PFCandidate cand;
+    //   FillCandTrack(cand, tracks[i]);
+    //   if( tkCaloMap.count(i)==0 ){
+    // 	unmatchedTks.push_back(cand);
+    //   } else { 
+    // 	// look for the first (highest-energy) unused matching cluster
+    // 	bool linked=false;
+    // 	for(int em_idx : tkCaloMap[i]){
+    // 	  if(!emUsed[em_idx]){
+    // 	    FillCandEMCalo(cand, tkCaloMap[i][0]);
+    // 	    caloMatchedTks.push_back(cand);
+    // 	    emUsed[ em_idx ] = true;
+    // 	    linked = true;
+    // 	    break;
+    // 	  }
+    // 	}
+    // 	if (!linked) unmatchedTks.push_back(cand);
+    //   }
+    // }
+    
+    // // ecal / hcal cluster arbitration
+    // std::vector<bool> hadUsed (hcalClusters.size(), false);
+    // for(int i=0; i<ecalClusters.size(); i++){
+    //   if( emHadCaloMap.count(i)==0 ){
+    // 	unmatchedTks.push_back(cand);
+    //   } else { 
+    // 	// look for the first (highest-energy) unused matching cluster
+    // 	bool linked=false;
+    // 	for(int em_idx : tkCaloMap[i]){
+    // 	  if(!emUsed[em_idx]){
+    // 	    FillCandEMCalo(cand, tkCaloMap[i][0]);
+    // 	    caloMatchedTks.push_back(cand);
+    // 	    emUsed[ em_idx ] = true;
+    // 	    linked = true;
+    // 	    break;
+    // 	  }
+    // 	}
+    // 	if (!linked) unmatchedTks.push_back(cand);
+    //   }
+    // }
+
+
+    // std::vector<ldmx::PFCandidate> unmatchedEMs;
+    // for(int i=0; i<ecalClusters.size(); i++){
+    //   if(emUsed[i]) continue;
+    //   ldmx::PFCandidate cand;
+    //   FillCandEMCalo(cand, ecalClusters[i]);
+    // }
+      
+	  
+    // 	}
+    // 	if( tkCaloMap[i].size()==1 ){
+    // 	if(!emUsed[ tkCaloMap[i][0] ]){
+    // 	  FillCandEMCalo(cand, tkCaloMap[i][0]);
+    // 	  caloMatchedTks.push_back(cand);
+    // 	  emUsed[ tkCaloMap[i][0] ] = true;
+    // 	}
+    //   }
+    //   } else if( tkCaloMap[i].size()==1 ){
+    // 	if(!emUsed[ tkCaloMap[i][0] ]){
+    // 	  FillCandEMCalo(cand, tkCaloMap[i][0]);
+    // 	  caloMatchedTks.push_back(cand);
+    // 	  emUsed[ tkCaloMap[i][0] ] = true;
+    // 	}
+    //   }
+    // }
+
+  } else {
+    // Single-particle builder
+    ldmx::PFCandidate pf;
+    int pid=0;
+    if (tracks.size()){
+      FillCandTrack(pf, tracks[0]);
+      pid += 1;
+    }
+    if (ecalClusters.size()){
+      FillCandEMCalo(pf, ecalClusters[0]);
+      pid += 2;
+    }
+    if (hcalClusters.size()){
+      FillCandHadCalo(pf, hcalClusters[0]);
+      pid += 4;
+    }
+    pf.setPID(pid);
+    pf.setEnergy( pf.getEcalEnergy() + pf.getHcalEnergy() );
+    pfCands.push_back(pf);  
   }
-  pf.setPID(pid);
-  pf.setEnergy( pf.getEcalEnergy() + pf.getHcalEnergy() );
-  pfCands.push_back(pf);
 
   // // Super-simple builder
   // for(int i=0; i<tracks.size(); i++){