diff --git a/sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.cc b/sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.cc
index eb9a3b810..66b668824 100644
--- a/sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.cc
+++ b/sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.cc
@@ -32,7 +32,6 @@ namespace blip {
       
 
     // Loop over cryostats
-    std::cout<<"NCryostats: "<<fGeom.Ncryostats()<<"\n";
     for(size_t cstat=0; cstat<fGeom.Ncryostats(); cstat++){
       auto const& cryoid = geo::CryostatID(cstat);
 
@@ -41,16 +40,14 @@ namespace blip {
         auto const& tpcid = geo::TPCID(cryoid,tpc);
 
         // Loop planes in TPC 'tpc'
-	auto const& plane0id = geo::PlaneID(cstat,tpc,0);
-	auto const& plane0geo = wireReadoutGeom->Get().Plane(plane0id);
+        auto const& plane0id = geo::PlaneID(cstat,tpc,0);
+        auto const& plane0geo = wireReadoutGeom->Get().Plane(plane0id);
         for(size_t pl=0; pl<wireReadoutGeom->Get().Nplanes(tpcid); pl++){
           auto const& planeid = geo::PlaneID(cstat,tpc,pl);
           auto const& planegeo = wireReadoutGeom->Get().Plane(planeid);
           kNumChannels += planegeo.Nwires();
             
           float offset = detProp.GetXTicksOffset(pl,tpc,cstat);
-          std::cout<<"CRYOSTAT "<<cstat<<" / TPC "<<tpc<<" / PLANE "<<pl<<":  "<<planegeo.Nwires()<<" wires\n";
-          std::cout<<"  XTicksOffset (from detProp): "<<offset<<"\n";
          
           kXTicksOffsets[cstat][tpc][pl] = 0;
 
@@ -66,12 +63,9 @@ namespace blip {
             auto const& tpcgeom   = cryostat.TPC(tpc);
             auto const xyz        = plane0geo.GetCenter();
             const double dir((tpcgeom.DriftSign() == geo::DriftSign::Negative) ? +1.0 : -1.0);
-	    float x_ticks_coefficient = kDriftVelocity*kTickPeriod;
+            float x_ticks_coefficient = kDriftVelocity*kTickPeriod;
             float goofy_offset = -xyz.X() / (dir * x_ticks_coefficient);
-	    std::cout<<"  After geometric correction: "<<offset - goofy_offset<<"\n";
-
             kXTicksOffsets[cstat][tpc][pl] = offset - goofy_offset;
-
 	  } else {
           
             // for the case of 2D wirecell workflow, the plane-to-plane
@@ -82,12 +76,7 @@ namespace blip {
             // we still want to correct for global trigger offset though:
             kXTicksOffsets[cstat][tpc][pl] = clockData.TriggerTime()/kTickPeriod;
 
-          }
-          
-          // additional ad-hoc corrections supplied by user
-          //kXTicksOffsets[cstat][tpc][pl] += fTimeOffset[pl];
-          std::cout << " offsetting plane " << pl << " by " << fTimeOffset[pl] << " ticks " << std::endl;
-          
+          }          
         }
       }
     }
@@ -164,18 +153,18 @@ namespace blip {
       for(int i=0; i<kNplanes; i++) {
         if( i == fCaloPlane ) continue;
         h_clust_overlap[iTPC][i]           = hdir.make<TH1D>(Form("t%i_p%i_clust_overlap",iTPC,i),   Form("TPC %i, Plane %i clusters;Overlap fraction",iTPC,i),101,0,1.01);
-	h_clust_dt[iTPC][i]                = hdir.make<TH1D>(Form("t%i_p%i_clust_dt",iTPC,i),        Form("TPC %i, Plane %i clusters;dT [ticks]",iTPC,i),200,-10,10);
+        h_clust_dt[iTPC][i]                = hdir.make<TH1D>(Form("t%i_p%i_clust_dt",iTPC,i),        Form("TPC %i, Plane %i clusters;dT [ticks]",iTPC,i),200,-10,10);
         h_clust_dtfrac[iTPC][i]            = hdir.make<TH1D>(Form("t%i_p%i_clust_dtfrac",iTPC,i),    Form("TPC %i, Plane %i clusters;Charge-weighted mean dT/RMS",iTPC,i),150,-1.5,1.5);
         
         h_clust_q[iTPC][i]     = hdir.make<TH2D>(Form("t%i_p%i_clust_charge",iTPC,i),  
 						 Form("Pre-cut, TPC %i;Plane %i cluster charge [#times10^{3} e-];Plane %i cluster charge [#times10^{3} e-]",iTPC,fCaloPlane,i),
 						 qbins,0,qmax,qbins,0,qmax);
-	h_clust_q[iTPC][i]->SetOption("colz");
+        h_clust_q[iTPC][i]->SetOption("colz");
         
         h_clust_q_cut[iTPC][i]     = hdir.make<TH2D>(Form("t%i_p%i_clust_charge_cut",iTPC,i),  
 						     Form("Post-cut, TPC %i;Plane %i cluster charge [#times10^{3} e-];Plane %i cluster charge [#times10^{3}]",iTPC,fCaloPlane,i),
 						     qbins,0,qmax,qbins,0,qmax);
-	h_clust_q_cut[iTPC][i]->SetOption("colz");
+        h_clust_q_cut[iTPC][i]->SetOption("colz");
       
         h_clust_score[iTPC][i]    = hdir.make<TH1D>(Form("t%i_p%i_clust_matchscore",iTPC,i),   Form("TPC %i, Plane %i clusters;Match score",iTPC,i),101,0,1.01);
        
@@ -189,7 +178,7 @@ namespace blip {
         h_clust_truematch_q[iTPC][i]     = hdir.make<TH2D>(Form("t%i_p%i_clust_truematch_charge",iTPC,i),  
           Form("Pre-cut, TPC %i;Plane %i cluster charge [#times10^{3} e-];Plane %i cluster charge [#times10^{3} e-]",iTPC,fCaloPlane,i),
           qbins,0,qmax,qbins,0,qmax);
-          h_clust_truematch_q[iTPC][i]->SetOption("colz");
+        h_clust_truematch_q[iTPC][i]->SetOption("colz");
         
         h_clust_truematch_score[iTPC][i]    = hdir.make<TH1D>(Form("t%i_p%i_clust_truematch_matchscore",iTPC,i),   Form("TPC %i, Plane %i clusters;Match score",iTPC,i),101,0,1.01);
 
@@ -239,6 +228,7 @@ namespace blip {
   void BlipRecoAlg::reconfigure( fhicl::ParameterSet const& pset ){
     
     fHitProducer        = pset.get<std::string>   ("HitProducer",       "gaushit");
+    fHitTruthMatcher    = pset.get<std::string>   ("HitTruthMatcher",   "blipgaushitTruthMatch");
     fTrkProducer        = pset.get<std::string>   ("TrkProducer",       "pandora");
     fGeantProducer      = pset.get<std::string>   ("GeantProducer",     "largeant");
     fSimDepProducer     = pset.get<std::string>   ("SimEDepProducer",   "ionandscint");
@@ -357,45 +347,39 @@ namespace blip {
     if (evt.getByLabel(fGeantProducer,pHandle))
       art::fill_ptr_vector(plist, pHandle);
  
-    // -- SimEnergyDeposits (usually dropped in reco
-    //art::Handle<std::vector<sim::EnergyDeposit> > sedHandle;
-    std::vector<sim::IDE > sedlist;
-    //if (evt.getByLabel(fSimDepProducer,sedHandle)){
-    //  art::fill_ptr_vector(sedlist, sedHandle);
-    // }
-    // -- SimChannels (usually dropped in reco)
+    // -- SimEnergyDeposits 
+    art::Handle<std::vector<sim::SimEnergyDeposit> > sedHandle;
+    std::vector<art::Ptr<sim::SimEnergyDeposit> > sedlist;
+    if (evt.getByLabel(fSimDepProducer,sedHandle)){
+      art::fill_ptr_vector(sedlist, sedHandle);
+     }
+    std::vector<sim::IDE > sIDElist;
+    // -- SimChannels 
     art::Handle<std::vector<sim::SimChannel> > simchanHandle;
     std::vector<art::Ptr<sim::SimChannel> > simchanlist;
     if (evt.getByLabel(fSimChanProducer,simchanHandle))
-      { 
+    { 
       art::fill_ptr_vector(simchanlist, simchanHandle);
-      //Loop over channels to get full sedlist
+      //Loop over channels to get full sIDElist
       for(int chIndex=0; chIndex<int(simchanlist.size()); chIndex++)
-	{
-	  std::vector<geo::WireID> wids    = wireReadoutGeom->Get().ChannelToWire( (*(simchanlist[chIndex])).Channel() ); //Not sure why this is a vector, but it should have len 1
-	  const geo::PlaneID&      planeID = wids[0].planeID();
-	  if(int(planeID.Plane) != fCaloPlane) continue; //only take calorimetry plane IDE values 
-	  std::vector< sim::IDE > TempChIDE = (*simchanlist[chIndex]).TrackIDsAndEnergies(0, 999999999);
-	  for(int ideIndex=0; ideIndex<int(TempChIDE.size()); ideIndex++)
 	    {
-	      //art::fill_ptr_vector(sedlist, simchanHandle.TrackIDsAndEnergies(0, 99999999));
-	      sedlist.push_back( TempChIDE[ideIndex] ); //may need to add a &
+      std::vector<geo::WireID> wids    = wireReadoutGeom->Get().ChannelToWire( (*(simchanlist[chIndex])).Channel() ); //Not sure why this is a vector, but it should have len 1
+      const geo::PlaneID&      planeID = wids[0].planeID();
+      if(int(planeID.Plane) != fCaloPlane) continue; //only take calorimetry plane IDE values 
+      unsigned int MaxTDCTick = UINT_MAX;
+      std::vector< sim::IDE > TempChIDE = (*simchanlist[chIndex]).TrackIDsAndEnergies(0, MaxTDCTick);
+      for(int ideIndex=0; ideIndex<int(TempChIDE.size()); ideIndex++)
+	        {
+	        sIDElist.push_back( TempChIDE[ideIndex] ); //may need to add a &
+	       }
 	    }
-	}
-      }
-
+    }
     // -- hits (from input module, usually track-masked subset of gaushit)
     art::Handle< std::vector<recob::Hit> > hitHandle;
     std::vector<art::Ptr<recob::Hit> > hitlist;
     if (evt.getByLabel(fHitProducer,hitHandle))
       art::fill_ptr_vector(hitlist, hitHandle);
 
-    // -- hits (from gaushit), these are used in truth-matching of hits
-    art::Handle< std::vector<recob::Hit> > hitHandleGH;
-    std::vector<art::Ptr<recob::Hit> > hitlistGH;
-    if (evt.getByLabel("gaushit",hitHandleGH))
-      art::fill_ptr_vector(hitlistGH, hitHandleGH);
-
     // -- tracks
     art::Handle< std::vector<recob::Track> > tracklistHandle;
     std::vector<art::Ptr<recob::Track> > tracklist;
@@ -404,8 +388,8 @@ namespace blip {
   
     // -- associations
     art::FindManyP<recob::Track> fmtrk(hitHandle,evt,fTrkProducer);
-    art::FindManyP<recob::Track> fmtrkGH(hitHandleGH,evt,fTrkProducer);
-    art::FindMany<simb::MCParticle,anab::BackTrackerHitMatchingData> fmhh(hitHandleGH,evt,"gaushitTruthMatch");
+    art::FindManyP<recob::Track> fmtrkGH(hitHandle,evt,fTrkProducer);
+    art::FindMany<simb::MCParticle, anab::BackTrackerHitMatchingData> fmhh(hitHandle,evt,fHitTruthMatcher);
     /*
     //====================================================
     // Update map of bad channels for this event
@@ -441,23 +425,7 @@ namespace blip {
     // use gaushitTruthMatch later on)
     //===============================================================
     std::map< int, int > map_gh;
-    // if input collection is already gaushit, this is trivial
-    if( fHitProducer == "gaushit" ) {
-      for(auto& h : hitlist ) map_gh[h.key()] = h.key(); 
-    // ... but if not, find the matching gaushit. There's no convenient
-    // hit ID, so we must loop through and compare channel/time (ugh)
-    } else {
-      std::map<int,std::vector<int>> map_chan_ghid;
-      for(auto& gh : hitlistGH ) map_chan_ghid[gh->Channel()].push_back(gh.key());
-      for(auto& h : hitlist ) {
-        for(auto& igh : map_chan_ghid[h->Channel()]){
-          if( hitlistGH[igh]->PeakTime() != h->PeakTime() ) continue;
-          map_gh[h.key()] = igh;
-          break;
-        }
-      }
-    }
-   
+    for(auto& h : hitlist ) map_gh[h.key()] = h.key(); 
     //=====================================================
     // Record PDG for every G4 Track ID
     //=====================================================
@@ -526,15 +494,23 @@ namespace blip {
     if( plist.size() ) {
       pinfo.resize(plist.size());
       for(size_t i = 0; i<plist.size(); i++){
-        BlipUtils::FillParticleInfo( *plist[i], pinfo[i], sedlist, fCaloPlane);
+        //use sim::EnergyDeposits by default. This is heavy and may be dropped
+        if(sedlist.size()>0)
+        {
+          BlipUtils::FillParticleInfo( *plist[i], pinfo[i], sedlist, fCaloPlane);
+        }
+        else //use sim::Channel -> IDE otherwise. This is usually kept but results in strange bugs.
+        {
+          BlipUtils::FillParticleInfo( *plist[i], pinfo[i], sIDElist, fCaloPlane);
+        }
         if( map_g4trkid_charge[pinfo[i].trackId] ) pinfo[i].numElectrons = (int)map_g4trkid_charge[pinfo[i].trackId];
+        else pinfo[i].numElectrons = pinfo[i].depElectrons; // JACOB ADDITION Jan22, 2026 for strange channel IDE results
         pinfo[i].index = i;
       }
       BlipUtils::MakeTrueBlips(pinfo, trueblips);
       BlipUtils::MergeTrueBlips(trueblips, fTrueBlipMergeDist);
     }
 
-
     //=======================================
     // Map track IDs to the index in the vector
     //=======================================
@@ -647,14 +623,8 @@ namespace blip {
       
 
       // find associated track
-      if( fHitProducer == "gaushit" && fmtrk.isValid() ) {
+      if( fmtrk.isValid() ) {
         if(fmtrk.at(i).size()) hitinfo[i].trkid = fmtrk.at(i)[0]->ID();
-      
-      // if the hit collection didn't have associations made
-      // to the tracks, try gaushit instead
-      } else if ( fmtrkGH.isValid() && map_gh.size() ) {
-        int gi = map_gh[i];
-        if (fmtrkGH.at(gi).size()) hitinfo[i].trkid= fmtrkGH.at(gi)[0]->ID(); 
       }
 
       // add to the map
@@ -664,8 +634,7 @@ namespace blip {
       if( hitinfo[i].trkid < 0 ) nhits_untracked++;
       //printf("  %lu   plane: %i,  wire: %i, time: %i\n",i,hitinfo[i].plane,hitinfo[i].wire,int(hitinfo[i].driftTime));
 
-    }//endloop over hits
-    
+    }//endloop over hits    
     //for(auto& a : tpc_plane_hitsMap ) {
       //for(auto& b : a.second ) 
         //std::cout<<"TPC "<<a.first<<", plane "<<b.first<<": "<<b.second.size()<<" hits\n";
@@ -690,6 +659,7 @@ namespace blip {
     
     
     // Basic track inclusion cut: exclude hits that were tracked
+    int Tracked=0;
     for(size_t i=0; i<hitlist.size(); i++){
       if( hitinfo[i].trkid < 0 ) continue;
       auto it = map_trkid_index.find(hitinfo[i].trkid);
@@ -698,9 +668,9 @@ namespace blip {
       if( tracklist[trkindex]->Length() > fMaxHitTrkLength ) {
         hitIsTracked[i] = true;
         hitIsGood[i] = false;
+        Tracked++;
       }
-    }
-        
+    }        
 
     // Filter based on hit properties. For hits that are a part of
     // multi-gaussian fits (multiplicity > 1), need to re-think this.
@@ -731,7 +701,6 @@ namespace blip {
     // Hit clustering
     // ---------------------------------------------------
     std::map<int,std::map<int,std::vector<int>>> tpc_planeclustsMap;
-   
     for(auto const& tpc_plane_hitsMap : cryo_tpc_plane_hitsMap ) {
     
     for(auto const& plane_hitsMap : tpc_plane_hitsMap.second ) {
@@ -892,7 +861,6 @@ namespace blip {
     
     float _matchQDiffLimit= (fMatchQDiffLimit <= 0 ) ? std::numeric_limits<float>::max() : fMatchQDiffLimit;
     float _matchMaxQRatio = (fMatchMaxQRatio  <= 0 ) ? std::numeric_limits<float>::max() : fMatchMaxQRatio;
-     
     for(auto& tpcMap : tpc_planeclustsMap ) { // loop on TPCs
       
       //std::cout
@@ -935,15 +903,15 @@ namespace blip {
               // Check that the two central wires intersect
               // *******************************************
               double y, z;
-	      geo::Point_t intsec_p;
-	      std::vector<geo::WireID> A_wireids = wireReadoutGeom->Get().ChannelToWire((unsigned int)hcA.CenterChan);
-	      std::vector<geo::WireID> B_wireids = wireReadoutGeom->Get().ChannelToWire((unsigned int)hcB.CenterChan);
+	            geo::Point_t intsec_p;
+	            std::vector<geo::WireID> A_wireids = wireReadoutGeom->Get().ChannelToWire((unsigned int)hcA.CenterChan);
+	            std::vector<geo::WireID> B_wireids = wireReadoutGeom->Get().ChannelToWire((unsigned int)hcB.CenterChan);
 
-	      if( !wireReadoutGeom->Get().WireIDsIntersect(A_wireids.at(0),B_wireids.at(0),intsec_p)) continue;
-	      // Save intersect location, so we don't have to
+	            if( !wireReadoutGeom->Get().WireIDsIntersect(A_wireids.at(0),B_wireids.at(0),intsec_p)) continue;
+	            // Save intersect location, so we don't have to
               // make another call to the Geometry service later
-	      y = intsec_p.Y();
-	      z = intsec_p.Z();
+	            y = intsec_p.Y();
+	            z = intsec_p.Z();
               TVector3 xloc(0,y,z);
               hcA.IntersectLocations[hcB.ID] = xloc;
               hcB.IntersectLocations[hcA.ID] = xloc;
@@ -1046,7 +1014,6 @@ namespace blip {
             for(auto& hc : hcGroup ) {
               hitclust[hc.ID].isMatched = true;
               for(auto hit : hitclust[hc.ID].HitIDs) hitinfo[hit].ismatch = true;
-             
               // Diagnostic plots for successful 3-plane matches
               //if( picky && hc.Plane != fCaloPlane ) {
                 //float q1 = (float)newBlip.clusters[fCaloPlane].Charge;
@@ -1096,18 +1063,26 @@ namespace blip {
             // if we made it this far, the blip is good!
             // associate this blip with the hits and clusters within it
             newBlip.ID = blips.size();
-            blips.push_back(newBlip);
             for(auto& hc : hcGroup ) {
               hitclust[hc.ID].BlipID = newBlip.ID;
               for( auto& h : hc.HitIDs ) hitinfo[h].blipid = newBlip.ID;
             }
+            //BLIPS HAVE A COPY OF HITCLUSTERS NOT A POINTER
+            //UPDATE THE HITCLUSTER VARS THAT HAVE CHANGED SINCE CONSTRUCTION
+            for(int iclust=0; iclust<int(sizeof(newBlip.clusters)/sizeof(newBlip.clusters[0])); iclust++)
+            {
+              if(newBlip.clusters[iclust].ID>-1){
+                newBlip.clusters[iclust].BlipID = newBlip.ID;
+                newBlip.clusters[iclust].isMatched = true;
+              }
+            }
+            blips.push_back(newBlip);
 
   
           }//endif ncands > 0
         }//endloop over caloplane ("Plane A") clusters
       }//endif calo plane has clusters
     }//endloop over TPCs
-
     // Re-index the clusters after removing unmatched
     if( !keepAllClusts ) {
       std::vector<blip::HitClust> hitclust_filt;
diff --git a/sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.h b/sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.h
index 59d125999..665fa30c2 100644
--- a/sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.h
+++ b/sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.h
@@ -129,6 +129,7 @@ namespace blip {
 
     // --- FCL configs ---
     std::string         fHitProducer;
+    std::string         fHitTruthMatcher;
     std::string         fTrkProducer;
     std::string         fGeantProducer;
     std::string         fSimDepProducer;
diff --git a/sbndcode/BlipRecoSBND/BlipAna_module.cc b/sbndcode/BlipRecoSBND/BlipAna_module.cc
index b9ca904e8..6f1b1c61d 100644
--- a/sbndcode/BlipRecoSBND/BlipAna_module.cc
+++ b/sbndcode/BlipRecoSBND/BlipAna_module.cc
@@ -456,10 +456,10 @@ class BlipAnaTreeDataStruct
       evtTree->Branch("hit_ismatch",hit_ismatch,"hit_ismatch[nhits]/I");
       evtTree->Branch("hit_trkid",hit_trkid,"hit_trkid[nhits]/I"); 
       if( saveTruthInfo ) {
-      evtTree->Branch("hit_g4trkid",hit_g4trkid,"hit_g4trkid[nhits]/I");
-      evtTree->Branch("hit_g4frac",hit_g4frac,"hit_g4frac[nhits]/F"); 
-      evtTree->Branch("hit_g4energy",hit_g4energy,"hit_g4energy[nhits]/F"); 
-      evtTree->Branch("hit_g4charge",hit_g4charge,"hit_g4charge[nhits]/F"); 
+        evtTree->Branch("hit_g4trkid",hit_g4trkid,"hit_g4trkid[nhits]/I");
+        evtTree->Branch("hit_g4frac",hit_g4frac,"hit_g4frac[nhits]/F"); 
+        evtTree->Branch("hit_g4energy",hit_g4energy,"hit_g4energy[nhits]/F"); 
+        evtTree->Branch("hit_g4charge",hit_g4charge,"hit_g4charge[nhits]/F"); 
       }
       evtTree->Branch("hit_clustid",hit_clustid,"hit_clustid[nhits]/I"); 
       evtTree->Branch("hit_blipid",hit_blipid,"hit_blipid[nhits]/I");
@@ -745,22 +745,22 @@ class BlipAna : public art::EDAnalyzer
     
     h_part_process    = dir_truth.make<TH1D>("part_process","MCParticle->Process()",5,0,5);
     auto xa = h_part_process->GetXaxis();
-      xa->SetBinLabel(1,"primary");
-      xa->SetBinLabel(2,"compt");
-      xa->SetBinLabel(3,"phot");
-      xa->SetBinLabel(4,"conv");
-      xa->SetBinLabel(5,"other");
+    xa->SetBinLabel(1,"primary");
+    xa->SetBinLabel(2,"compt");
+    xa->SetBinLabel(3,"phot");
+    xa->SetBinLabel(4,"conv");
+    xa->SetBinLabel(5,"other");
     
     h_nblips_tm    = dir_truth.make<TH1D>("nblips_tm","Truth-matched 3D blips per event",blipBins,0,blipMax);
     h_blip_qcomp   = dir_truth.make<TH1D>("blip_qcomp","Fraction of true charge (at anode) reconstructed into 3D blips",202,0,1.01);
     h_blip_reszy   = dir_truth.make<TH2D>("blip_res_zy","Blip position resolution;Z_{reco} - Z_{true} [cm];Y_{reco} - Y_{true} [cm]",150,-15,15,150,-15,15);
-      h_blip_reszy ->SetOption("colz");
+    h_blip_reszy ->SetOption("colz");
     h_blip_resx    = dir_truth.make<TH1D>("blip_res_x","Blip position resolution;X_{reco} - X_{true} [cm]",200,-10,10);
     h_blip_resE   = dir_truth.make<TH1D>("blip_res_E","Blip energy resolution;(E_{reco} - E_{true})/E_{true} [cm]",200,-1.,1.);
     h_blip_E_vs_resE    = dir_truth.make<TH2D>("blip_res_energy","Energy resolution of 3D blips;Energy [MeV];#deltaE/E_{true}",100,0,5,200,-1.,1.);
-        h_blip_E_vs_resE ->SetOption("colz");
+    h_blip_E_vs_resE ->SetOption("colz");
     h_clust_qres_vs_q       = dir_truth.make<TH2D>("qres_vs_q","Clusters on collection plane;True charge deposited [ #times 10^{3} e- ];Reco resolution",160,0,80,200,-1,1);
-      h_clust_qres_vs_q     ->SetOption("colz");
+    h_clust_qres_vs_q     ->SetOption("colz");
     h_clust_qres_anode      = dir_truth.make<TH1D>("qres_anode","Reco charge vs true charge collected;( reco-true ) / true;Area-normalized entries",200,-1.,1.);
     h_clust_qres_dep        = dir_truth.make<TH1D>("qres_dep","Reco charge vs true charge deposited;( reco-true ) / true;Area-normalized entries",200,-1.,1.);
     h_qratio_vs_time_sim  = dir_truth.make<TH2D>("qratio_vs_time_sim",";Drift time [#mus]; Q_{anode} / Q_{dep}",44,100,2300, 1000,0.50,1.50);
@@ -801,10 +801,10 @@ class BlipAna : public art::EDAnalyzer
       h_hitqres[i]      = dir_truth.make<TH1D>(Form("pl%i_hit_q_res",i),     Form("Plane %i hits;charge resolution: (reco-true)/true",i),    400,-1,1);
       h_hitqres_scatter[i] = dir_truth.make<TH2D>( Form("pl%i_hit_qres_scatter",i), 
         Form("Plane %i;true hit charge [#times 10^{3} e-];Reconstructed hit charge [#times 10^{3} e-]",i),qBins,0,qMax/1e3,qBins,0,qMax/1e3);
-        h_hitqres_scatter[i]  ->SetOption("colz");
+      h_hitqres_scatter[i]  ->SetOption("colz");
       h_hitqres_vs_q[i] = dir_truth.make<TH2D>( Form("pl%i_hit_qres_vs_q",i),
         Form("Plane %i;true hit charge [#times 10^{3} e-];hit charge resolution: (reco-true)/true",i),qBins,0,qMax/1e3, 400,-2,2);
-        h_hitqres_vs_q[i]     ->SetOption("colz");
+      h_hitqres_vs_q[i]     ->SetOption("colz");
     
       h_chargecomp[i] = dir_truth.make<TH1D>(Form("pl%i_hit_charge_completeness",i),Form("charge completness, plane %i",i),101,0,1.01);
       h_hitpur[i]     = dir_truth.make<TH1D>(Form("pl%i_hit_purity",i),Form("hit purity, plane %i",i),101,0,1.01);
@@ -1326,7 +1326,6 @@ void BlipAna::analyze(const art::Event& evt)
       }
     
     }
-      
     if( fDebugMode ) PrintClusterInfo(clust);
     
   }//endloop over 2D hit clusters
@@ -1386,7 +1385,7 @@ void BlipAna::analyze(const art::Event& evt)
       nblips_picky++;
       fNum3DBlipsPicky++;
       h_blip_charge_picky ->Fill(blp.clusters[fCaloPlane].Charge);
-       h_blip_zy_picky     ->Fill(blp.Position.Z(), blp.Position.Y());
+      h_blip_zy_picky     ->Fill(blp.Position.Z(), blp.Position.Y());
       h_blip_charge_YU_picky->Fill( 0.001*blp.clusters[2].Charge, 0.001*blp.clusters[0].Charge );
       h_blip_charge_YV_picky->Fill( 0.001*blp.clusters[2].Charge, 0.001*blp.clusters[1].Charge );
       h_blip_charge_UV_picky->Fill( 0.001*blp.clusters[0].Charge, 0.001*blp.clusters[1].Charge );
@@ -1489,27 +1488,27 @@ void BlipAna::endJob(){
   //printf("                 picky frac   : %5.3f\n",     fNum3DBlipsPicky/float(fNum3DBlips));
   
   if(fIsMC){
-  printf("  MC-matched blips per evt    : %.3f\n",       fNum3DBlipsTrue/nEvents);
-  printf("  MC blip purity              : %.3f\n",       fNum3DBlipsTrue/float(fNum3DBlips));
-  printf("  MC blip purity, 3 planes    : %.3f\n",      fNum3DBlipsTrue3P/float(fNum3DBlips3Plane));
-  if( h_blip_qcomp->GetMean() > 0 ) 
-  printf("  Charge completeness, total  : %.4f +/- %.4f\n", h_blip_qcomp->GetMean(), h_blip_qcomp->GetStdDev()/sqrt(fNumEvents));
-  //printf("                       < 2MeV : %.4f +/- %.4f\n", h_blip_qcomp_2MeV->GetMean(), h_blip_qcomp_2MeV->GetStdDev()/sqrt(fNumEvents));
-  //printf("  Blip purity                 : %.4f\n",       h_blip_pur->GetMean());
+    printf("  MC-matched blips per evt    : %.3f\n",       fNum3DBlipsTrue/nEvents);
+    printf("  MC blip purity              : %.3f\n",       fNum3DBlipsTrue/float(fNum3DBlips));
+    printf("  MC blip purity, 3 planes    : %.3f\n",      fNum3DBlipsTrue3P/float(fNum3DBlips3Plane));
+    if( h_blip_qcomp->GetMean() > 0 ) 
+    printf("  Charge completeness, total  : %.4f +/- %.4f\n", h_blip_qcomp->GetMean(), h_blip_qcomp->GetStdDev()/sqrt(fNumEvents));
+    //printf("                       < 2MeV : %.4f +/- %.4f\n", h_blip_qcomp_2MeV->GetMean(), h_blip_qcomp_2MeV->GetStdDev()/sqrt(fNumEvents));
+    //printf("  Blip purity                 : %.4f\n",       h_blip_pur->GetMean());
   }
   printf("  Mean blip charge            : %.0f e-\n",      h_blip_charge->GetMean());
   printf("\n");
   for(size_t i=0; i<kNplanes; i++){
-  printf("  Plane %lu -------------------------\n",i);
-  printf("   * total hits/evt           : %.2f\n",fNumHits[i]/(float)fNumEvents);
-  printf("   * untracked hits/evt       : %.2f (%.2f plane-matched)\n",fNumHitsUntracked[i]/(float)fNumEvents, fNumHitsMatched[i]/(float)fNumEvents);
-  //printf("   * plane-matched hits/evt   : %.2f\n",fNumHitsMatched[i]/(float)fNumEvents);
-  if(fIsMC) {
-  printf("   * true-matched hits/evt    : %.2f (%.2f plane-matched)\n",fNumHitsTrue[i]/(float)fNumEvents, fNumHitsMatchedTrue[i]/(float)fNumEvents);
-  if( h_chargecomp[i]->GetMean() > 0 ) 
-  printf("   * charge completeness      : %.4f\n",h_chargecomp[i]->GetMean());
-  printf("   * hit purity               : %.4f\n",h_hitpur[i]->GetMean());
-  }
+    printf("  Plane %lu -------------------------\n",i);
+    printf("   * total hits/evt           : %.2f\n",fNumHits[i]/(float)fNumEvents);
+    printf("   * untracked hits/evt       : %.2f (%.2f plane-matched)\n",fNumHitsUntracked[i]/(float)fNumEvents, fNumHitsMatched[i]/(float)fNumEvents);
+    //printf("   * plane-matched hits/evt   : %.2f\n",fNumHitsMatched[i]/(float)fNumEvents);
+    if(fIsMC) {
+      printf("   * true-matched hits/evt    : %.2f (%.2f plane-matched)\n",fNumHitsTrue[i]/(float)fNumEvents, fNumHitsMatchedTrue[i]/(float)fNumEvents);
+      if( h_chargecomp[i]->GetMean() > 0 ) 
+      printf("   * charge completeness      : %.4f\n",h_chargecomp[i]->GetMean());
+      printf("   * hit purity               : %.4f\n",h_hitpur[i]->GetMean());
+    }
   } 
   printf("\n***********************************************\n");
 
@@ -1582,6 +1581,9 @@ void BlipAna::PrintClusterInfo(const blip::HitClust& hc){
     hc.EdepID,
     hc.isMatched
   );
+  printf("G4 IDs contib \n");
+  for(int g4ID : hc.G4IDs) printf(" %i ", g4ID);
+  printf("\n");
 }
 
 void BlipAna::PrintBlipInfo(const blip::Blip& bl){
diff --git a/sbndcode/BlipRecoSBND/Utils/BlipUtils.cc b/sbndcode/BlipRecoSBND/Utils/BlipUtils.cc
index e1b8b1b48..8c1dd53af 100644
--- a/sbndcode/BlipRecoSBND/Utils/BlipUtils.cc
+++ b/sbndcode/BlipRecoSBND/Utils/BlipUtils.cc
@@ -32,9 +32,8 @@ namespace BlipUtils {
   //===========================================================================
   // Provided a MCParticle, calculate everything we'll need for later calculations
   // and save into ParticleInfo object
-  void FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo, SEDVec_t& sedvec, int caloPlane){
-    
-    // Get important info and do conversions
+void FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo){
+  // Get important info and do conversions
     pinfo.particle    = part;
     pinfo.trackId     = part.TrackId();
     pinfo.isPrimary   = (int)(part.Process() == "primary");
@@ -50,17 +49,28 @@ namespace BlipUtils {
     pinfo.time        = /*ns ->mus*/1e-3 * part.T();
     pinfo.endtime     = /*ns ->mus*/1e-3 * part.EndT();
     pinfo.numTrajPts  = part.NumberTrajectoryPoints();
-
     // Pathlength (in AV) and start/end point
     pinfo.pathLength  = PathLength( part, pinfo.startPoint, pinfo.endPoint);
-
     // Central position of trajectory
     pinfo.position    = 0.5*(pinfo.startPoint+pinfo.endPoint);
-
     // Energy/charge deposited by this particle, found using SimEnergyDeposits 
     pinfo.depEnergy     = 0;
     pinfo.depElectrons  = 0;
+    return;
+}
+void FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo, SEDVec_t& sedvec, int caloPlane){
+    FillParticleInfo( part, pinfo);
     for(auto& sed : sedvec ) {
+      if( -1*sed->TrackID() == part.TrackId() || sed->TrackID() == part.TrackId() ) {
+        pinfo.depEnergy     += sed->Energy();
+        pinfo.depElectrons  += sed->NumElectrons();
+      }
+    }
+    return;
+  }
+  void FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo, SIDEVec_t& sIDEvec, int caloPlane){
+    FillParticleInfo( part, pinfo);
+    for(auto& sed : sIDEvec ) {
       if( -1*sed.trackID == part.TrackId() || sed.trackID == part.TrackId() ) {
         pinfo.depEnergy     += sed.energy;
         pinfo.depElectrons  += sed.numElectrons;
@@ -99,12 +109,14 @@ namespace BlipUtils {
 
       // We want to loop through any contiguous electrons (produced
       // with process "eIoni") and add the energy they deposit into this blip.
-      if( part.NumberDaughters() ) {
+      if( part.NumberDaughters() ) { //particles have daughters but they must all be neutron, gamma, or one of the special processes?
         for(size_t j=0; j<pinfo.size(); j++){
           simb::MCParticle& p = pinfo[j].particle;
           std::string pr = p.Process();
-          if( p.PdgCode() != 2112 && (pr == "eIoni" || pr == "muIoni" || pr == "hIoni") ){
-            if( IsAncestorOf(p.TrackId(),part.TrackId(),true) ) GrowTrueBlip(pinfo[j],tb);
+          if( p.PdgCode() != 2112 && p.PdgCode() != 22 && (pr == "eIoni" || pr == "muIoni" || pr == "hIoni") ){ //neutron and photons leave track
+            if( IsAncestorOf(p.TrackId(),part.TrackId(),true,true) ){
+              GrowTrueBlip(pinfo[j],tb);
+            }
           }
         }
       }
@@ -151,7 +163,7 @@ namespace BlipUtils {
     // .. otherwise, check that the new particle
     // creation time is comparable to existing blip.
     // then calculate new energy-weighted position.
-    } else if ( fabs(tblip.Time-pinfo.time) < 3 ) {
+    } else if ( fabs(tblip.Time-pinfo.time) < 3) {
       float totE = tblip.Energy + pinfo.depEnergy;
       float w1 = tblip.Energy/totE;
       float w2 = pinfo.depEnergy/totE;
@@ -383,7 +395,7 @@ namespace BlipUtils {
       // use view with the maximal wire extent to calculate transverse (YZ) length
       if( hcs[i].NWires > newblip.MaxWireSpan ) {
         newblip.MaxWireSpan = hcs[i].NWires;
-	newblip.dYZ         = hcs[i].NWires * wirepitch;
+	      newblip.dYZ         = hcs[i].NWires * wirepitch;
       }
   
       for(size_t j=i+1; j<hcs.size(); j++){
@@ -397,9 +409,8 @@ namespace BlipUtils {
           intsec_p.SetY(hcs[i].IntersectLocations.find(hcs[j].ID)->second.Y());
           intsec_p.SetZ(hcs[i].IntersectLocations.find(hcs[j].ID)->second.Z());
         } else {
-	  std::vector<geo::WireID> i_wireids = wireReadoutGeom->Get().ChannelToWire((unsigned int)hcs[i].CenterChan);
-	  std::vector<geo::WireID> j_wireids = wireReadoutGeom->Get().ChannelToWire((unsigned int)hcs[j].CenterChan);
-
+	        std::vector<geo::WireID> i_wireids = wireReadoutGeom->Get().ChannelToWire((unsigned int)hcs[i].CenterChan);
+	        std::vector<geo::WireID> j_wireids = wireReadoutGeom->Get().ChannelToWire((unsigned int)hcs[j].CenterChan);
           match3d = wireReadoutGeom->Get().WireIDsIntersect(i_wireids.at(0), j_wireids.at(0), intsec_p);
         }
 
@@ -474,7 +485,7 @@ namespace BlipUtils {
   //====================================================================
   // Function to determine if a particle descended from another particle.
   // Allows option to break lineage at photons for contiguous parentage.
-  bool IsAncestorOf(int particleID, int ancestorID, bool breakAtPhots = false){
+  bool IsAncestorOf(int particleID, int ancestorID, bool breakAtPhots = false, bool breakAtNeutrons = false){
     art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
     const sim::ParticleList& plist = pi_serv->ParticleList();
     if( particleID == ancestorID  )       return true;
@@ -488,6 +499,7 @@ namespace BlipUtils {
       simb::MCParticle pM = pi_serv->TrackIdToParticle(p.Mother());
       if      ( pM.TrackId() == ancestorID )                      { return true;  }
       else if ( breakAtPhots == true && pM.PdgCode() == 22 )      { return false; }
+      else if ( breakAtNeutrons && pM.PdgCode() == 2112 )           { return false; }
       else if ( pM.Process() == "primary" || pM.TrackId() == 1 )  { return false; }
       else if ( pM.Mother() == 0 )                                { return false; }
       else    { particleID = pM.TrackId(); }              
diff --git a/sbndcode/BlipRecoSBND/Utils/BlipUtils.h b/sbndcode/BlipRecoSBND/Utils/BlipUtils.h
index 0584c9792..59ed79735 100644
--- a/sbndcode/BlipRecoSBND/Utils/BlipUtils.h
+++ b/sbndcode/BlipRecoSBND/Utils/BlipUtils.h
@@ -37,8 +37,8 @@
 #include "sbndcode/BlipRecoSBND/Utils/DataTypes.h"
 #include "TH1D.h"
 
-
-typedef std::vector<sim::IDE> SEDVec_t;
+typedef std::vector<art::Ptr<sim::SimEnergyDeposit>> SEDVec_t;
+typedef std::vector<sim::IDE> SIDEVec_t;
 
 geo::View_t kViews[3]={geo::kU, geo::kV, geo::kW};
 
@@ -48,7 +48,9 @@ namespace BlipUtils {
   // Functions related to blip reconstruction
   //###################################################
   //void      InitializeDetProps();
-  void      FillParticleInfo(simb::MCParticle const&, blip::ParticleInfo&, SEDVec_t&, int plane=2);
+  void      FillParticleInfo(simb::MCParticle const&, blip::ParticleInfo&, SEDVec_t&, int plane);
+  void      FillParticleInfo(simb::MCParticle const&, blip::ParticleInfo&, SIDEVec_t&, int plane);
+  void      FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo);
   //void      CalcPartDrift(blip::ParticleInfo&, int);
   //void      CalcTotalDep(float&,int&,float&, SEDVec_t&);
   void      MakeTrueBlips(std::vector<blip::ParticleInfo>&, std::vector<blip::TrueBlip>&); 
@@ -78,7 +80,7 @@ namespace BlipUtils {
   //bool    G4IdToMCTruth( int const, art::Ptr<simb::MCTruth>&);
   double  PathLength(const simb::MCParticle&, TVector3&, TVector3&);
   double  PathLength(const simb::MCParticle&);
-  bool    IsAncestorOf(int, int, bool);
+  bool    IsAncestorOf(int, int, bool, bool);
   double  DistToBoundary(const recob::Track::Point_t&);
   double  DistToLine(TVector3&, TVector3&, TVector3&);
   double  DistToLine2D(TVector2&, TVector2&, TVector2&);
diff --git a/sbndcode/BlipRecoSBND/blipreco_configs.fcl b/sbndcode/BlipRecoSBND/blipreco_configs.fcl
index 1d4158a21..01329758a 100644
--- a/sbndcode/BlipRecoSBND/blipreco_configs.fcl
+++ b/sbndcode/BlipRecoSBND/blipreco_configs.fcl
@@ -6,9 +6,10 @@ BEGIN_PROLOG
 sbnd_blipalg:
 {
   HitProducer:          "specialblipgaushit"       #// input recob::Hits to use for blip reconstruction
+  HitTruthMatcher:      "blipgaushitTruthMatch"    #// Truth info corresponding to the above hit collection
   TrkProducer:          "blipPandoraTrackCopy"       #// input recob::Tracks to use for blip reconstruction
   GeantProducer:        "largeant"      #// input sim::MCParticles (getting true particle info)
-  SimEDepProducer:      "ionandscint"   #// input sim::SimEnergyDeposits (getting energy/electrons deposited)
+  SimEDepProducer:      "ionandscint:priorSCE:G4"   #// input sim::SimEnergyDeposits (getting energy/electrons deposited)
   SimChanProducer:      "simtpc2d:simpleSC:DetSim"      #// label for sim::SimChannels (getting drifted charge; optional)
   MaxHitTrkLength:      5.0             #// hits in trks > this length will be vetoed [cm]
   DoHitFiltering:       false           #// filter hits based on amp, width, RMS