diff --git a/CMakeLists.txt b/CMakeLists.txt
index 06c9db9f..723321bd 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -29,6 +29,7 @@ find_package(larreco REQUIRED)
 find_package(nugen REQUIRED)
 find_package(nutools REQUIRED)
 find_package(Eigen3 REQUIRED)
+find_package( TorchScatter QUIET EXPORT)
 
 cet_cmake_module_directories(Modules BINARY)
 
diff --git a/ubreco/BlipReco/Alg/BlipRecoAlg.cc b/ubreco/BlipReco/Alg/BlipRecoAlg.cc
index 485829ba..fa87f34c 100644
--- a/ubreco/BlipReco/Alg/BlipRecoAlg.cc
+++ b/ubreco/BlipReco/Alg/BlipRecoAlg.cc
@@ -1,4 +1,5 @@
 #include "ubreco/BlipReco/Alg/BlipRecoAlg.h"
+#include "lardataobj/AnalysisBase/Calorimetry.h"
 #include <fstream>
 
 namespace blip {
@@ -19,7 +20,7 @@ namespace blip {
     kTickPeriod           = clockData.TPCClock().TickPeriod();
     kNominalRecombFactor  = ModBoxRecomb(fCalodEdx,kNominalEfield);
     kWion                 = 1000./util::kGeVToElectrons;
-    
+
     // initialize channel list
     fBadChanMask       .resize(8256,false);
     fBadChanMaskPerEvt = fBadChanMask;
@@ -45,7 +46,7 @@ namespace blip {
     int NBadChansFromFile     = std::count(fBadChanMask.begin(),fBadChanMask.end(),true);
     
     EvtBadChanCount = 0;
-
+   
     printf("******************************************\n");
     printf("Initializing BlipRecoAlg...\n");
     printf("  - Efield: %.4f kV/cm\n",detProp.Efield());
@@ -147,18 +148,22 @@ namespace blip {
   //###########################################################
   void BlipRecoAlg::reconfigure( fhicl::ParameterSet const& pset ){
     
+    fDebugMode  = pset.get<bool> ("DebugMode",false);
+
     // initialize MC flags (will be determined automatically
     // later after the first event is processed)
     isMC        = false;
     isMCOverlay = false;
   
     fHitProducer        = pset.get<art::InputTag> ("HitProducer","");
+    fHitProducerOG      = pset.get<art::InputTag> ("HitProducerOG","");
     fHitProducerData    = pset.get<art::InputTag> ("HitProducerData",    "gaushit::DataRecoStage1Test");
     fHitProducerOverlay = pset.get<art::InputTag> ("HitProducerOverlay", "gaushit::OverlayStage1a");
-    fHitProducerOG      = fHitProducerData;
-    if( fHitProducer=="" ) fHitProducer = fHitProducerData;
-    fHitTruthMatch      = pset.get<art::InputTag> ("HitTruthMatch",     "gaushitTruthMatch::OverlayRecoStage1b");
+    if( fHitProducerOG=="" ) fHitProducerOG = fHitProducerData;
+    if( fHitProducer=="" )   fHitProducer   = fHitProducerData;
+    fHitTruthMatch      = pset.get<art::InputTag> ("HitTruthMatch",     ""); //"gaushitTruthMatch::OverlayRecoStage1b");
     fTrkProducer        = pset.get<std::string>   ("TrkProducer",       "pandoraInit");
+    fTrkProducerCal     = pset.get<std::string>   ("TrkProducerCal",    "pandoracaliInit");
     fGeantProducer      = pset.get<std::string>   ("GeantProducer",     "largeant");
     fSimDepProducer     = pset.get<std::string>   ("SimEDepProducer",   "ionization");
     fSimChanProducer    = pset.get<std::string>   ("SimChanProducer",   "driftWC:simpleSC");
@@ -193,8 +198,8 @@ namespace blip {
     
     fMinMatchedPlanes   = pset.get<int>           ("MinMatchedPlanes",    2);
     fPickyBlips         = pset.get<bool>          ("PickyBlips",          false);
-    fApplyTrkCylinderCut= pset.get<bool>          ("ApplyTrkCylinderCut", false);
-    fCylinderRadius     = pset.get<float>         ("CylinderRadius",      15);
+    //fApplyTrkCylinderCut= pset.get<bool>          ("ApplyTrkCylinderCut", false);
+    //fCylinderRadius     = pset.get<float>         ("CylinderRadius",      15);
     fIgnoreDataTrks     = pset.get<bool>          ("IgnoreDataTrks",      false);
 
     fCaloPlane          = pset.get<int>           ("CaloPlane",           2);
@@ -210,7 +215,7 @@ namespace blip {
     fBadChanProducer    = pset.get<std::string>   ("BadChanProducer",     "nfspl1:badchannels");
     fBadChanFile        = pset.get<std::string>   ("BadChanFile",         "");
     fMinDeadWireGap     = pset.get<int>           ("MinDeadWireGap",      1);
-    
+   
     //fKeepAllClusts[0] = pset.get<bool>          ("KeepAllClustersInd", false);
     //fKeepAllClusts[1] = pset.get<bool>          ("KeepAllClustersInd", false);
     //fKeepAllClusts[2] = pset.get<bool>          ("KeepAllClustersCol", true);
@@ -276,12 +281,15 @@ namespace blip {
     //=====================================================
     // Record PDG for every G4 Track ID
     //=====================================================
-    //std::map<int,int> map_g4trkid_pdg;
+    map_g4trkid_index.clear();
     map_g4trkid_pdg.clear();
     map_g4trkid_chan.clear();
     map_g4trkid_chan_energy.clear();
     map_g4trkid_chan_charge.clear();
-    for(size_t i = 0; i<plist.size(); i++) map_g4trkid_pdg[plist[i]->TrackId()] = plist[i]->PdgCode();
+    for(size_t i = 0; i<plist.size(); i++) {
+      map_g4trkid_pdg[plist[i]->TrackId()] = plist[i]->PdgCode();
+      map_g4trkid_index[plist[i]->TrackId()] = i;
+    }
     //std::map<int, std::set<int>>         map_g4trkid_chan;
     //std::map<int, std::map<int,double> > map_g4trkid_chan_energy;
     //std::map<int, std::map<int,double> > map_g4trkid_chan_charge;
@@ -365,10 +373,7 @@ namespace blip {
       }
     }
 
-    
-
-
-  }
+  }// RunBlipTruth
   
   
 
@@ -380,10 +385,12 @@ namespace blip {
   // this function is run.
   //###########################################################
   void BlipRecoAlg::RunBlipReco( const art::Event& evt ) {
-  
-    //std::cout<<"\n"
-    //<<"=========== BlipRecoAlg =========================\n"
-    //<<"Event "<<evt.id().event()<<" / run "<<evt.id().run()<<"\n";
+    
+    if( fDebugMode ) {
+      std::cout<<"\n"
+      <<"=========== BlipRecoAlg =========================\n"
+      <<"Event "<<evt.id().event()<<" / run "<<evt.id().run()<<"\n";
+    }
     
     //=======================================
     // Extract truth info if applicable
@@ -403,8 +410,14 @@ namespace blip {
     hitclust.clear();
     hitinfo.clear();
     EvtBadChanCount = 0;
-    
   
+    map_blip_hitids.clear();
+    map_blip_trkID.clear();
+    map_blip_trkIDfrac.clear();
+    map_blip_primaryPDG.clear();
+    map_blip_primaryG4ID.clear();
+    map_blip_ncategory.clear();
+
     //=======================================
     // Get data products for this event
     //========================================
@@ -416,6 +429,9 @@ namespace blip {
     auto const& tpcCalib_provider   = art::ServiceHandle<lariov::TPCEnergyCalibService>()->GetProvider();
     auto const& chanFilt            = art::ServiceHandle<lariov::ChannelStatusService>()->GetProvider();
     
+    // lifetime (ms)
+    kLifetime = lifetime_provider.Lifetime();
+    
     //====================================================
     // Update map of bad channels for this event
     //====================================================
@@ -457,7 +473,8 @@ namespace blip {
     // -- associations
     art::FindManyP<recob::Track> fmtrk(hitHandle,evt,fTrkProducer);
     art::FindManyP<recob::Track> fmtrkOG(hitHandleOG,evt,fTrkProducer);
-    
+    art::FindManyP<anab::Calorimetry> fmcal(tracklistHandle,evt,fTrkProducerCal);
+
     //===============================================================
     // Map of each hit to its gaushit index (needed if the provided
     // hit collection is some filtered subset of gaushit, in order to
@@ -485,7 +502,7 @@ namespace blip {
     //=======================================
     // Map track IDs to the index in the vector
     //=======================================
-    //std::cout<<"Looping over tracks...\n";
+    if( fDebugMode ) std::cout<<"Looping over "<<tracklist.size()<<" tracks...\n";
     //std::map<size_t,size_t> map_trkid_isMC;
     map_trkid_isMC.clear();
     map_trkid_index.clear();
@@ -508,7 +525,7 @@ namespace blip {
     std::map<int,std::vector<int>> planehitsMap;
     int nhits_untracked = 0;
 
-    //std::cout<<"Looping over the hits... "<<hitlist.size()<<"\n";
+    if( fDebugMode ) std::cout<<"Looping over "<<hitlist.size()<<" hits...\n";
     for(size_t i=0; i<hitlist.size(); i++){
       auto const& thisHit = hitlist[i];
       int   chan    = thisHit->Channel();
@@ -558,9 +575,46 @@ namespace blip {
     //========================================
     if( pinfo.size() ) {
       
+      // If specific association tag label was not provided, check which ones are 
+      // in the event and identify which is most compatible with the hit collection. 
+      // This allows (for example) the wiremod workflow to not require a separate FCL
+      // configuration file specifying a different association label.
+      //
+      // If a valid association is found, fHitTruthMatch is updated, so this check
+      // is done only once. This assumes the data files being analyzed all have consistent
+      // data product labels.
+      //
+      if( fHitTruthMatch.label() == "" ) {
+
+        int best_match_count = 0;
+
+        // Find all MCParticle <-> Hit association products in the event
+        auto tags = evt.getInputTags<art::Assns<simb::MCParticle,recob::Hit,anab::BackTrackerHitMatchingData>>();
+        for (auto &tag : tags) {
+          if( tag.label() != "gaushitTruthMatch") continue;
+          art::FindMany<simb::MCParticle,anab::BackTrackerHitMatchingData> fmhh(hitHandleOG,evt,tag);
+
+          // is it valid?
+          if( !fmhh.isValid() ) continue;
+
+          // how many matches are there?
+          int nmatches = 0;
+          for(size_t i=0; i<hitlist.size(); i++) nmatches+=fmhh.at(map_gh[i]).size();
+          if( nmatches <= best_match_count ) continue;
+
+          // if we are good so far, update the tag 
+          best_match_count = nmatches;
+          fHitTruthMatch = tag;
+        }
+
+        if (fDebugMode) std::cout<<"*** Found valid MCP <--> hit assn:  "<<fHitTruthMatch.label()<<"::"<<fHitTruthMatch.process()<<" ***\n";      
+
+      }
+      if (fDebugMode) std::cout<<"Doing hit truth-matching with assn: "<<fHitTruthMatch.label()<<"::"<<fHitTruthMatch.process()<<"\n";      
+
       // Get backtracker associations
       art::FindMany<simb::MCParticle,anab::BackTrackerHitMatchingData> fmhh(hitHandleOG,evt,fHitTruthMatch);
-    
+
       // Loop the hits
       for(size_t i=0; i<hitlist.size(); i++){
         
@@ -570,40 +624,44 @@ namespace blip {
         // the truth-matching metadata is stored in the event
         //--------------------------------------------------
         int igh = map_gh[i];
-        if( fmhh.at(igh).size() ) {
-          std::vector<simb::MCParticle const*> pvec;
-          std::vector<anab::BackTrackerHitMatchingData const*> btvec;
-          fmhh.get(igh,pvec,btvec);
-          hitinfo[i].g4energy = 0;
-          hitinfo[i].g4charge = 0;
-          float maxQ = -9;
-          for(size_t j=0; j<pvec.size(); j++){
-            hitinfo[i].g4energy += btvec.at(j)->energy;
-            hitinfo[i].g4charge += btvec.at(j)->numElectrons;
-            if( btvec.at(j)->numElectrons <= maxQ ) continue;
-            maxQ = btvec.at(j)->numElectrons;
-            hitinfo[i].g4trkid  = pvec.at(j)->TrackId();
-            hitinfo[i].g4pdg    = pvec.at(j)->PdgCode();
-            hitinfo[i].g4frac   = btvec.at(j)->ideNFraction;
-          }
+        if( fmhh.isValid() ) {
           
-          // ###      uB behavior as of Nov 2022              ###
-          // WireCell's detsim implements its gain "fudge factor" 
-          // by scaling the SimChannel electrons instead of the gain
-          // response. So we need to correct for this effect to get 
-          // accurate count of 'true' electrons collected on channel.
-          if( fSimGainFactor > 0 ) hitinfo[i].g4charge /= fSimGainFactor;
-         
-          /* 
-          // Some old code to measure wire-by-wire thresholds
-          if( map_g4trkid_chan_energy[hitinfo[i].g4trkid][chan] > 0 ) {
-            double trueEnergyDep = map_g4trkid_chan_energy[hitinfo[i].g4trkid][chan];
-            h_recoWireEff_num->Fill(trueEnergyDep);}
-          if( map_g4trkid_chan_charge[hitinfo[i].g4trkid][chan] > 0 ) {
-            double trueChargeDep = map_g4trkid_chan_charge[hitinfo[i].g4trkid][chan];
-            h_recoWireEffQ_num->Fill(trueChargeDep);}
-          */
-        }
+          if( fmhh.at(igh).size() ) {
+            std::vector<simb::MCParticle const*> pvec;
+            std::vector<anab::BackTrackerHitMatchingData const*> btvec;
+            fmhh.get(igh,pvec,btvec);
+            hitinfo[i].g4energy = 0;
+            hitinfo[i].g4charge = 0;
+            float maxQ = -9;
+            for(size_t j=0; j<pvec.size(); j++){
+              hitinfo[i].g4energy += btvec.at(j)->energy;
+              hitinfo[i].g4charge += btvec.at(j)->numElectrons;
+              if( btvec.at(j)->numElectrons <= maxQ ) continue;
+              maxQ = btvec.at(j)->numElectrons;
+              hitinfo[i].g4trkid  = pvec.at(j)->TrackId();
+              hitinfo[i].g4pdg    = pvec.at(j)->PdgCode();
+              hitinfo[i].g4frac   = btvec.at(j)->ideNFraction;
+            }
+            
+            // ###      uB behavior as of Nov 2022              ###
+            // WireCell's detsim implements its gain "fudge factor" 
+            // by scaling the SimChannel electrons instead of the gain
+            // response. So we need to correct for this effect to get 
+            // accurate count of 'true' electrons collected on channel.
+            if( fSimGainFactor > 0 ) hitinfo[i].g4charge /= fSimGainFactor;
+           
+            /* 
+            // Some old code to measure wire-by-wire thresholds
+            if( map_g4trkid_chan_energy[hitinfo[i].g4trkid][chan] > 0 ) {
+              double trueEnergyDep = map_g4trkid_chan_energy[hitinfo[i].g4trkid][chan];
+              h_recoWireEff_num->Fill(trueEnergyDep);}
+            if( map_g4trkid_chan_charge[hitinfo[i].g4trkid][chan] > 0 ) {
+              double trueChargeDep = map_g4trkid_chan_charge[hitinfo[i].g4trkid][chan];
+              h_recoWireEffQ_num->Fill(trueChargeDep);}
+            */
+          }
+
+        }//fmhh isValid
         
         // IF this hit was (a) matched to a track, and (b) matched to a truth
         // energy deposit, then keep the tally
@@ -856,7 +914,7 @@ namespace blip {
       
       }
     }
-
+    if( fDebugMode ) std::cout<<"Created "<<hitclust.size()<<" hit clusters...\n";
 
     // =============================================================================
     // Plane matching and 3D blip formation
@@ -1041,46 +1099,82 @@ namespace blip {
 
             
             // ----------------------------------------
-            // apply cylinder cut 
-            for(auto& trk : tracklist ){
+            // Loop over all >5cm tracks in the event, find point
+            // of closest approach
+            //for(auto& trk : tracklist ){
+            for (size_t itrk = 0; itrk<tracklist.size(); itrk++){
+              auto &trk = tracklist[itrk];
               if( trk->Length() < fMaxHitTrkLength ) continue;
               if( fIgnoreDataTrks && !map_trkid_isMC[trk->ID()] ) continue;
               auto& a = trk->Vertex();
               auto& b = trk->End();
               TVector3 p1(a.X(), a.Y(), a.Z() );
               TVector3 p2(b.X(), b.Y(), b.Z() );
-              // TO-DO: if this track starts or ends at a TPC boundary, 
-              // we should extend p1 or p2 to outside the AV to avoid blind spots
-              TVector3 bp = newBlip.Position;
-              float dToLine = BlipUtils::DistToLine(p1,p2,bp);
-              float d = dToLine;
-              if( dToLine < 0 ) d = std::min( (bp-p1).Mag(), (bp-p2).Mag() );
-              if( d > 0 ) {
+              float L = (p2-p1).Mag();
+              
+              float d = -9;
+
+              // First determine if the blip is closer to an interior point
+              // of the track, or to one of the endpoints
+              float d1 = (newBlip.Position-p1).Mag();
+              float d2 = (newBlip.Position-p2).Mag();
+              if( fabs( pow(d1,2) - pow(d2,2) ) > pow(L,2) ){
+                d = std::min(d1,d2);
+              
+              // if closer to an interior point, check track spacepoints
+              // and look for closest approach
+              } else if (fmcal.isValid()) {
+                float closestApproach = 9999;
+                //float maxDeviation = 0;
+                std::vector<art::Ptr<anab::Calorimetry> > caloObjs = fmcal.at(itrk);
+                for(auto &caloObj : caloObjs ) {
+                  //std::cout<<"  found "<<caloObj->XYZ().size()<<" 3D points for a track length of "<<trk->Length()<<" ("<<caloObj->XYZ().size()/trk->Length()<<")\n";
+                  for(auto &xyz : caloObj->XYZ() ) {
+                    TVector3 sp(xyz.X(),xyz.Y(),xyz.Z());
+                    float dist = (newBlip.Position-sp).Mag();
+                    if( dist < closestApproach) closestApproach = dist;
+                    //float deviation = BlipUtils::DistToLine(p1,p2,sp);
+                    //if( deviation > maxDeviation ) maxDeviation = deviation;
+                  }
+                }
                 
-                // update closest trkdist
-                if( newBlip.ProxTrkDist < 0 || d < newBlip.ProxTrkDist ) {
+                d = closestApproach;
+              
+
+              // if no spacepoints available, just project a line between start/end
+              // of the track and calculate the point of closest approach algebraically
+              } else {
+                d = BlipUtils::DistToLine(p1,p2,newBlip.Position);
+              }
+              
+              if (d<0) continue;
+
+              // update closest trkdist
+              if( newBlip.ProxTrkDist < 0 || d < newBlip.ProxTrkDist ) {
                   newBlip.ProxTrkDist = d;
                   newBlip.ProxTrkID = trk->ID();
-                }
+              }
+                
                 // need to do some math to figure out if this is in
                 // the 45 degreee "cone" relative to the start/end
-                if( dToLine > 0 && !newBlip.inCylinder && d < fCylinderRadius ) {
-                  float angle1 = asin( d / (p1-bp).Mag() ) * 180./3.14159;
-                  float angle2 = asin( d / (p2-bp).Mag() ) * 180./3.14159;
-                  if( angle1 < 45. && angle2 < 45. ) newBlip.inCylinder = true;
-                }
-              }
+                //if( dToLine > 0 && !newBlip.inCylinder && d < fCylinderRadius ) {
+                //  float angle1 = asin( d / (p1-bp).Mag() ) * 180./3.14159;
+                //  float angle2 = asin( d / (p2-bp).Mag() ) * 180./3.14159;
+                //  if( angle1 < 45. && angle2 < 45. ) newBlip.inCylinder = true;
+                //}
+              
+            
             }//endloop over trks
-           
-            if( fApplyTrkCylinderCut && newBlip.inCylinder ) continue;
+            
+            //if( fApplyTrkCylinderCut && newBlip.inCylinder ) continue;
             
             // In the case that a cluster appeared to be touching a track in 
             // one of the 3 views, but after 3D evaluation is found to be positioned
             // far away from that track, we must update the "TouchTrk" status.
             if( (newBlip.ProxTrkID != newBlip.TouchTrkID )
-            //    || (newBlip.ProxTrkDist > 1.5*newBlip.dX ) 
-            //    || (newBlip.ProxTrkDist > 1.5*newBlip.dYZ ) ) 
-            ) {
+                || (newBlip.ProxTrkDist > 1.5*newBlip.dX ) 
+                || (newBlip.ProxTrkDist > 1.5*newBlip.dYZ ) ) 
+            {
               newBlip.TouchTrkID = -9;
             }
 
@@ -1091,7 +1185,10 @@ namespace blip {
             blips.push_back(newBlip);
             for(auto& hc : hcGroup ) {
               hitclust[hc.ID].BlipID = newBlip.ID;
-              for( auto& h : hc.HitIDs ) hitinfo[h].blipid = newBlip.ID;
+              for( auto h : hc.HitIDs ) {
+                hitinfo[h].blipid = newBlip.ID;
+                map_blip_hitids[newBlip.ID].push_back(h);
+              }
             }
 
   
@@ -1099,6 +1196,8 @@ namespace blip {
         }//endloop over caloplane ("Plane A") clusters
       }//endif calo plane has clusters
     }//endloop over TPCs
+    
+    if( fDebugMode ) std::cout<<"Created "<<blips.size()<<" 3D blips...\n";
 
     // Re-index the clusters after removing unmatched
     //if( !keepAllClusts ) {
@@ -1121,7 +1220,6 @@ namespace blip {
     for(size_t i=0; i<hitlist.size(); i++){
       if (hitinfo[i].trkid >= 0 ) continue;
       h_chan_nhits->Fill(wireReadout.PlaneWireToChannel(geo::WireID(0, 0, hitinfo[i].plane, hitinfo[i].wire)));
-      
       int clustid = hitinfo[i].clustid;
       if( clustid >= 0 ) {
         if( hitclust[clustid].NWires > 1 ) continue;
@@ -1132,7 +1230,159 @@ namespace blip {
       //h_chan_nclusts->Fill(wireReadout.PlaneWireToChannel(hitinfo[i].plane,hitinfo[i].wire));
     }
 
+
+
+
+    //*************************************************************************
+    // Additional blip-level calculations
+    //*************************************************************************
+    for(size_t i=0; i<blips.size(); i++){
+      auto& blip = blips[i];
+     
+
+      //=====================================================
+      // Determine if blip shares hits with any tracks
+      //=====================================================
+      map_blip_trkID[blip.ID] = -9;
+      map_blip_trkIDfrac[blip.ID] = -9;
+      int nhits = map_blip_hitids[blip.ID].size();
+      int nbtids = 0;
+      if( nhits==0 ) continue;
+      std::map<int,int> map_btid_count;
+      for(auto bh : map_blip_hitids[blip.ID] ) {
+        if( hitinfo[bh].trkid < 0 ) continue;
+        nbtids++;
+        map_btid_count[hitinfo[bh].trkid]++;
+      }
+      auto max_val_it = std::max_element(std::begin(map_btid_count), std::end(map_btid_count),
+        [](const std::pair<int, int>& p1, const std::pair<int, int>& p2) { return p1.second < p2.second;});
+      if (nbtids>0 && max_val_it != std::end(map_btid_count)) {
+        map_blip_trkID[blip.ID] = max_val_it->first;
+        map_blip_trkIDfrac[blip.ID] = float(max_val_it->second) / nhits;
+      }
+     
+
+      // ================================================
+      // Save the true blip into the object;
+      // each cluster must match to the same energy dep
+      // ================================================
+      std::set<int> set_edepids;
+      for(auto& hc : blip.clusters ) {
+        if( !hc.isValid ) continue; 
+        if( hc.EdepID < 0 ) break;
+        set_edepids.insert( hc.EdepID );
+      }
+     
+      // initialize some dummy values
+      map_blip_primaryPDG[blip.ID]  = 0;
+      map_blip_primaryG4ID[blip.ID] = -99;
+      map_blip_ncategory[blip.ID]   = -9;
+      blip.truth.LeadG4PDG    = 0;
+      blip.truth.LeadG4ID     = -99;
+      blip.truth.LeadG4Index  = -99;
+
+      if( set_edepids.size() == 1 ){
+        
+        // save the true blip into the object
+        blip.truth = trueblips[*set_edepids.begin()];
+
+        // trace ancestry of this blip
+        auto part = BlipUtils::ReturnMCParticle(blip.truth.LeadG4ID);
+        int blipPDG = part.PdgCode();
+        auto blipProc = part.Process();
+        auto ancestors = BlipUtils::ReturnAllAncestors(blip.truth.LeadG4ID);
+        //std::cout<<"Blip "<<blip.ID<<" (PDG "<<part.PdgCode()<<", "<<part.Process()<<") has "<<ancestors.size()<<" ancestors: \n";
+        for(auto& p : ancestors) {
+          //std::cout<<"   "<<p.TrackId()<<"    PDG "<<p.PdgCode()<<"   process? "<<p.Process()<<"\n";
+          if( p.Process() == "primary" ) { 
+            map_blip_primaryPDG[blip.ID] = p.PdgCode();
+            map_blip_primaryG4ID[blip.ID] = p.TrackId();
+          }
+        }
+       
+        ///========================================
+        // Categories for neutron-based studies
+        //=========================================
+        //  -9 = no truth match (data/overlay)
+        //  0  = truth-matched, but not falling in category
+        //  1  = primary (n,1p)
+        //  2  = primary (n,Np)
+        //  3  = secondary (n,1p)
+        //  4  = secondary (n,Np)
+        //  5  = primary (n,gamma)
+        //  6  = secondary (n,gamma)
+        //  7  = ncapture gamma
+        //  8  = gamma from mu capture
+        if( fDebugMode ) std::cout<<"Determining neutron ancestors... blip PDG is "<<blipPDG<<"   with process "<<blipProc<<"\n";
+        if( ancestors.size() > 0 ) {
+          map_blip_ncategory[blip.ID] = 0;
+          
+          for(auto p : ancestors ) {
+            double x1 = p.Vx();
+            double y1 = p.Vy();
+            double z1 = p.Vz();
+            double x2 = p.EndX();
+            double y2 = p.EndY();
+            double z2 = p.EndZ();
+            if(fDebugMode) std::cout<<"      - G4TrkID: "<<p.TrackId()<<"     PDG: "<<p.PdgCode()<<"    Process: "<<p.Process()
+                                     <<", start: ("<<x1<<", "<<y1<<", "<<z1<<")   end: ("<<x2<<", "<<y2<<", "<<z2<<")\n";
+          }
+
+          //----------------------------------------------
+          // inelastic (n,p) - categories 1-4
+          bool primaryneut = false;
+          if( blipPDG == 2212 && blipProc == "neutronInelastic"){
+            if( ancestors[0].Process() == "primary" ) primaryneut = true; 
+            // figure out how many daughter protons this neutron made 
+            int np = 0;
+            int mothern = ancestors[0].TrackId();
+            for(size_t i=map_g4trkid_index[mothern]+1; i<pinfo.size(); i++){
+              if( pinfo[i].particle.Mother() == mothern && pinfo[i].particle.PdgCode() == 2212 ) np++;
+            }
+            //std::cout<<"Blip "<<blip.ID<<" came from immediate ancestor "<<ancestors[0].PdgCode()<<", that made "<<np<<" protons    "<<primaryneut<<"\n";
+            if      ( primaryneut && np == 1 )  map_blip_ncategory[blip.ID] = 1;
+            else if ( primaryneut && np > 1  )  map_blip_ncategory[blip.ID] = 2;
+            else if ( !primaryneut && np == 1 ) map_blip_ncategory[blip.ID] = 3; 
+            else if ( !primaryneut && np > 1 )  map_blip_ncategory[blip.ID] = 4;
+          } 
+          
+          //----------------------------------------------
+          // electron-based blips from gammas
+          else if ( abs(blipPDG) == 11 && blipProc != "primary" && ancestors.size() >= 2 ) {
+            
+            // if the immediate ancestor is a photon...
+            if( ancestors[0].PdgCode() == 22 ) {
+              
+              // Loop back over ancestors
+              for(size_t k=0; k<ancestors.size(); k++){
+                auto& a = ancestors[k];
+                if( a.Process() == "nCapture" ) { 
+                  map_blip_ncategory[blip.ID] = 7; 
+                  break;
+                }
+                if( a.PdgCode() == 22 && a.Process() == "muMinusCaptureAtRest" ) {
+                  map_blip_ncategory[blip.ID] = 8; 
+                  break;
+                }
+                if( a.PdgCode() == 22 && a.Process() == "neutronInelastic" ) { 
+                  map_blip_ncategory[blip.ID] = (ancestors[k+1].Process() == "primary") ? 5 : 6; 
+                  break;
+                }
+              }
+            }
+          }//endif blip type
+          
+          if(fDebugMode) std::cout<<"      category assigned: "<< map_blip_ncategory[blip.ID] <<"\n";
+
+        }//end blip ancestry
+        //std::cout<<"BLIP CATEGORY = "<<map_blip_ncategory[blip.ID]
+        //<<", primary ancestor PDG "<<map_blip_primaryPDG[blip.ID]<<" -------------------------\n";
+
+      }//endif blip has true edep 
     
+    }//endloop over blips
+
+
     //*************************************************************************
     // Loop over the vector of blips and perform calorimetry calculation.
     // Here we fill:
@@ -1166,10 +1416,10 @@ namespace blip {
       // --- Lifetime correction ---
       // Note: Without knowing real T0 of a blip, this correction is meaningless.
       //       Units of 'ms', not microseconds, hence the 1E-3 conversion factor.
-      if( fLifetimeCorr && blip.Time>0 ) {
-        float t = blip.Time*1e-3;
-        float tau = lifetime_provider.Lifetime();
-        blip.ChargeCorr = std::max(0.,(double)blip.Charge) * exp( t / tau ); 
+      if( fLifetimeCorr && blip.Time>0 && kLifetime > 0 ) {
+        //float t = blip.Time*1e-3;
+        //float tau = lifetime_provider.Lifetime();
+        blip.ChargeCorr = std::max(0.,(double)blip.Charge) * exp( blip.Time*1e-3 / kLifetime ); 
       }
       
       // --- SCE corrections ---
@@ -1216,22 +1466,11 @@ namespace blip {
       // METHOD 2 (TODO)
       //std::cout<<"Calculating ESTAR energy dep...  "<<depEl<<", "<<Efield<<"\n";
       //blips[i].EnergyESTAR = ESTAR->Interpolate(depEl, Efield); 
-      
-      // ================================================
-      // Save the true blip into the object;
-      // each cluster must match to the same energy dep
-      // ================================================
-      std::set<int> set_edepids;
-      for(auto& hc : blip.clusters ) {
-        if( !hc.isValid ) continue; 
-        if( hc.EdepID < 0 ) break;
-        set_edepids.insert( hc.EdepID );
-      }
-      if( set_edepids.size() == 1 )
-        blip.truth = trueblips[*set_edepids.begin()];
-      
     
     }//endloop over blip vector
+    
+
+
 
   }//End main blip reco function
  
@@ -1274,8 +1513,8 @@ namespace blip {
     printf("  Charge ratio maximum      : %.1f\n",        fMatchMaxQRatio);    
     printf("  Minimum match score       : %.2f\n",        fMatchMinScore);
     printf("  Ignoring data tracks?     : %i\n",          fIgnoreDataTrks);
-    printf("  Track-cylinder radius     : %.1f cm\n",     fCylinderRadius);
-    printf("  Applying cylinder cut?    : %i\n",          fApplyTrkCylinderCut);
+    //printf("  Track-cylinder radius     : %.1f cm\n",     fCylinderRadius);
+    //printf("  Applying cylinder cut?    : %i\n",          fApplyTrkCylinderCut);
     
     /*
     printf("  Min cluster overlap       : ");
diff --git a/ubreco/BlipReco/Alg/BlipRecoAlg.h b/ubreco/BlipReco/Alg/BlipRecoAlg.h
index 2e974a83..69681eac 100644
--- a/ubreco/BlipReco/Alg/BlipRecoAlg.h
+++ b/ubreco/BlipReco/Alg/BlipRecoAlg.h
@@ -63,7 +63,6 @@
 #include <math.h>
 #include <limits>
 
-
 namespace blip {
 
   //--------------------------------------------
@@ -108,7 +107,8 @@ namespace blip {
     //TH1D*   h_recoWireEff_num;
     //TH1D*   h_recoWireEffQ_denom;
     //TH1D*   h_recoWireEffQ_num;
-    
+   
+    std::map<int,int>                     map_g4trkid_index;
     std::map<int,int>                     map_g4trkid_pdg;
     std::map<int, std::set<int>>          map_g4trkid_chan;
     std::map<int, std::map<int,double> >  map_g4trkid_chan_energy;
@@ -120,10 +120,23 @@ namespace blip {
     float   kNominalEfield;
     float   kDriftVelocity;
     float   kTickPeriod;
-    int     kNumChannels;
+    //int     kNumChannels;
     float   kLifetime;
 
     art::InputTag     fHitProducer;
+    
+    // Map of blip hit IDs    
+    std::map<int,std::vector<int>> map_blip_hitids;
+    
+    // TrackID for every blip (if there's a pandora track that shares hits)
+    std::map<int,int>   map_blip_trkID;
+    std::map<int,float> map_blip_trkIDfrac;
+  
+    // Primary particle that blip descends from
+    std::map<int,int>   map_blip_primaryPDG;
+    std::map<int,int>   map_blip_primaryG4ID;
+    std::map<int,int>   map_blip_ncategory;
+    
 
    private:
       
@@ -133,11 +146,12 @@ namespace blip {
     art::InputTag       fHitProducerOG;
     art::InputTag       fHitTruthMatch;
     std::string         fTrkProducer;
+    std::string         fTrkProducerCal;
     std::string         fGeantProducer;
     std::string         fSimDepProducer;
     std::string         fSimChanProducer;
     float               fSimGainFactor;
-    //bool                fDebugMode;
+    bool                fDebugMode;
     float               fTrueBlipMergeDist;
     bool                fDoHitFiltering;
     float               fMaxHitTrkLength;
@@ -177,6 +191,8 @@ namespace blip {
     
     //bool                keepAllClusts;
     //bool                fKeepAllClusts[kNplanes];
+   
+    //blip::BlipRNNAlg*   fRNNAlg;
 
     // --- Calorimetry configs ---
     int                 fCaloPlane;
@@ -187,7 +203,6 @@ namespace blip {
     float               fModBoxA;
     float               fModBoxB;
     
- 
 
     // --- Histograms ---
     //TH1D*   h_chanstatus;
diff --git a/ubreco/BlipReco/BlipAnaReader_module.cc b/ubreco/BlipReco/BlipAnaReader_module.cc
index 49b8c0ed..a49ccf05 100644
--- a/ubreco/BlipReco/BlipAnaReader_module.cc
+++ b/ubreco/BlipReco/BlipAnaReader_module.cc
@@ -36,6 +36,9 @@
 #include "ubevt/Database/UbooneElectronLifetimeService.h"
 #include "larevt/SpaceChargeServices/SpaceChargeService.h"
 #include "ubreco/BlipReco/Alg/BlipRecoAlg.h"
+#include "ubreco/BlipReco/Utils/NuSelectionToolBase.h"
+#include "ubreco/BlipReco/Utils/NuSelectionSCECorrections.h"
+#include "ubreco/BlipReco/Utils/NuSelectionTrackShowerScoreFuncs.h"
 
 // C++ includes
 #include <cstring>
@@ -101,6 +104,15 @@ class BlipAnaReaderTreeDataStruct
   int           badchans;             // #bad chans according to wirecell
   int           longtrks;             // tracks > 5 cm
   
+  // --- MCTruth neutrino info ---
+  int   mctruth_nu_pdg;         // Neutrino PDG (if present, otherwise = 0)
+  int   mctruth_nu_ccnc;        // CC (0) or NC (1)
+  int   mctruth_nu_mode;        // interaction mode from Genie
+  float mctruth_nu_vtx_x;       // vertex X 
+  float mctruth_nu_vtx_y;       // vertex Y
+  float mctruth_nu_vtx_z;       // vertex Z
+  float mctruth_nu_KE;          // kinetic energy 
+  
   // --- 3D Blip information ---
   int   nblips;                       // number of blips in event
   int   blip_id[kMaxBlips];           // blip ID / index
@@ -127,6 +139,17 @@ class BlipAnaReaderTreeDataStruct
   int   blip_touchtrkid[kMaxBlips];   // track ID of touched track
   int   blip_clustid[kNplanes][kMaxBlips];     // cluster ID per plane
   
+  // --- Reconstructed neutrino slice information (Pandora) --
+  bool      nu_isNeutrino;            // neutrino slice identified by Pandora
+  vfloat_t  nu_nuscore;               // neutrino score
+  int       nu_pfp_pdg;               // PDG particle best matching with reco slice
+  float     nu_reco_vtx_x;            // reconstructed vertex X [cm]
+  float     nu_reco_vtx_y;            // reconstructed vertex Y [cm]
+  float     nu_reco_vtx_z;            // reconstructed vertex Z [cm]
+  vint_t    nu_trk_id;        // trackIDs for tracks in this PFP
+  vfloat_t  nu_trk_score;     // track scores for tracks in this PFP
+  vfloat_t  nu_shwr_score;    // shower scores in this PFP
+  
   // === Function for resetting data ===
   void Clear(){ 
     event                 = -999; // --- event-wide info ---
@@ -136,7 +159,23 @@ class BlipAnaReaderTreeDataStruct
     badchans              = -99;
     longtrks              = -99;
     timestamp             = -999;
-    //timestamp_hr          = -999;
+    
+    mctruth_nu_pdg        = 0;
+    mctruth_nu_ccnc       = -9;
+    mctruth_nu_mode       = -9;
+    mctruth_nu_vtx_x      = -999;
+    mctruth_nu_vtx_y      = -999;
+    mctruth_nu_vtx_z      = -999;
+    mctruth_nu_KE         = -999;
+    nu_isNeutrino         = false;
+    nu_nuscore            .clear();
+    nu_pfp_pdg            = -999;
+    nu_reco_vtx_x         = -999;
+    nu_reco_vtx_y         = -999;
+    nu_reco_vtx_z         = -999;
+    nu_trk_score          .clear();
+    nu_trk_id             .clear();
+    nu_shwr_score         .clear();  
    
     nblips                    = 0;
     FillWith(blip_id,         -9);
@@ -164,6 +203,8 @@ class BlipAnaReaderTreeDataStruct
   void MakeTree(){
     
     art::ServiceHandle<art::TFileService> tfs;
+    auto vf = "std::vector<float>";
+    auto vi = "std::vector<int>";
    
     evtTree = tfs->make<TTree>(treeName.c_str(),"analysis tree");
     evtTree->Branch("event",&event,"event/I");
@@ -174,6 +215,23 @@ class BlipAnaReaderTreeDataStruct
     evtTree->Branch("lifetime",&lifetime,"lifetime/F");
     evtTree->Branch("badchans",&badchans,"badchans/I");
     evtTree->Branch("longtrks",&longtrks,"longtrks/I");
+    evtTree->Branch("mctruth_nu_pdg",&mctruth_nu_pdg,"mctruth_nu_pdg/I");
+    evtTree->Branch("mctruth_nu_ccnc",&mctruth_nu_ccnc,"mctruth_nu_ccnc/I");
+    evtTree->Branch("mctruth_nu_mode",&mctruth_nu_mode,"mctruth_nu_mode/I");
+    evtTree->Branch("mctruth_nu_vtx_x",&mctruth_nu_vtx_x,"mctruth_nu_vtx_x/F");
+    evtTree->Branch("mctruth_nu_vtx_y",&mctruth_nu_vtx_y,"mctruth_nu_vtx_y/F");
+    evtTree->Branch("mctruth_nu_vtx_z",&mctruth_nu_vtx_z,"mctruth_nu_vtx_z/F");
+    evtTree->Branch("mctruth_nu_KE",&mctruth_nu_KE,"mctruth_nu_KE/F");
+    evtTree->Branch("nu_isNeutrino",&nu_isNeutrino,"nu_isNeutrino/O");
+    //evtTree->Branch("nu_nuscore",&nu_nuscore,"nu_nuscore/F");
+    evtTree->Branch("nu_nuscore", vf, &nu_nuscore);
+    evtTree->Branch("nu_pfp_pdg",&nu_pfp_pdg,"nu_pfp_pdg/I");
+    evtTree->Branch("nu_reco_vtx_x",&nu_reco_vtx_x,"nu_reco_vtx_x/F");
+    evtTree->Branch("nu_reco_vtx_y",&nu_reco_vtx_y,"nu_reco_vtx_y/F");
+    evtTree->Branch("nu_reco_vtx_z",&nu_reco_vtx_z,"nu_reco_vtx_z/F");
+    evtTree->Branch("nu_trk_id", vi, &nu_trk_id);
+    evtTree->Branch("nu_trk_score", vf, &nu_trk_score);
+    evtTree->Branch("nu_shwr_score", vf, &nu_shwr_score);
     evtTree->Branch("nblips",&nblips,"nblips/I");
     evtTree->Branch("blip_nplanes",blip_nplanes,"blip_nplanes[nblips]/I");
     evtTree->Branch("blip_x",blip_x,"blip_x[nblips]/F");
@@ -229,6 +287,27 @@ class BlipAnaReader : public art::EDAnalyzer
   int   fNum3DBlips         = 0;
   int   fNum3DBlips3Plane   = 0;
   int   fNum3DBlipsTrue     = 0;
+  
+  // --- Neutrino selection tools
+  using ProxyPfpColl_t = selection::ProxyPfpColl_t;
+  using ProxyPfpElem_t = selection::ProxyPfpElem_t;
+  art::InputTag fPFPproducer;
+  art::InputTag fCLSproducer; // cluster associated to PFP
+  art::InputTag fSLCproducer; // slice associated to PFP
+  art::InputTag fHITproducer; // hit associated to cluster
+  art::InputTag fSHRproducer; // shower associated to PFP
+  art::InputTag fVTXproducer; // vertex associated to PFP
+  art::InputTag fPCAproducer; // PCAxis associated to PFP
+  art::InputTag fMCTproducer;
+  art::InputTag fTRKproducer;
+
+  //void    BuildPFPMap(const ProxyPfpColl_t&);
+  //void    AddDaughters(const ProxyPfpElem_t &pfp_pxy,
+  //                  const ProxyPfpColl_t &pfp_pxy_col,
+  //                  std::vector<ProxyPfpElem_t> &slice_v);
+
+  // a map linking the PFP Self() attribute used for hierarchy building to the PFP index in the event record
+  std::map<unsigned int, unsigned int> _pfpmap;
 
   // --- Histograms ---
   TH1D*   h_nblips;
@@ -290,9 +369,20 @@ BlipAnaReader::BlipAnaReader(fhicl::ParameterSet const& pset) :
   fData = new BlipAnaReaderTreeDataStruct();
   fData ->Clear();
   fData ->MakeTree();
-
+  
   // initialize histograms
   InitializeHistograms();
+  
+  fPFPproducer = pset.get<art::InputTag>("PFPproducer","pandora");
+  fSHRproducer = pset.get<art::InputTag>("SHRproducer","shrreco3d");
+  fHITproducer = pset.get<art::InputTag>("HITproducer","pandora");
+  fVTXproducer = pset.get<art::InputTag>("VTXproducer","pandora");
+  fPCAproducer = pset.get<art::InputTag>("PCAproducer","pandora");
+  fCLSproducer = pset.get<art::InputTag>("CLSproducer","pandora");
+  fSLCproducer = pset.get<art::InputTag>("SLCproducer","pandora");
+  fMCTproducer = pset.get<art::InputTag>("MCTproducer","generator");
+  fTRKproducer = pset.get<art::InputTag>("TRKproducer","pandora");
+
     
 }
 BlipAnaReader::~BlipAnaReader(){}
@@ -332,6 +422,140 @@ void BlipAnaReader::analyze(const art::Event& evt)
   <<"Event "<<evt.id().event()<<" / run "<<evt.id().run()<<"; total: "<<fNumEvents<<"\n";
   }
   
+  //===================================
+  // Check neutrinos in MCTruth
+  // (NuanceOffset found through simb::kNuanceOffset)
+  //===================================
+  art::Handle< std::vector<simb::MCTruth> > truthHandle;
+  std::vector<art::Ptr<simb::MCTruth> > truthlist;
+  if (evt.getByLabel("generator",truthHandle))
+    art::fill_ptr_vector(truthlist, truthHandle); 
+  if( truthlist.size() > 0 ) {
+    auto& nu    = truthlist[0]->GetNeutrino();
+    auto& part  = truthlist[0]->GetParticle(0);
+    auto PDG    = fabs(part.PdgCode());
+    if ( (PDG == 12) || (PDG == 14) ) {
+      fData->mctruth_nu_pdg   = part.PdgCode();
+      fData->mctruth_nu_ccnc  = nu.CCNC();
+      fData->mctruth_nu_mode  = nu.Mode();
+      fData->mctruth_nu_vtx_x = part.EndPosition()[0];
+      fData->mctruth_nu_vtx_y = part.EndPosition()[1];
+      fData->mctruth_nu_vtx_z = part.EndPosition()[2];
+      fData->mctruth_nu_KE    =  /*GeV->MeV*/1e3 * (part.E()-part.Mass());
+    }
+  }
+  
+  
+  //=======================================================
+  // Check if the event contains neutrino information
+  // (much of this copied from ubana/searchingfornues)
+  //=======================================================
+  // -- PFPs
+  art::Handle< std::vector<recob::PFParticle> > pfpHandle;
+  std::vector<art::Ptr<recob::PFParticle> > pfplist;
+  if (evt.getByLabel(fPFPproducer,pfpHandle))
+    art::fill_ptr_vector(pfplist, pfpHandle);
+
+  // -- associated tracks/vertex
+  if( pfplist.size() ) {
+    // grab PFParticles in event
+    ProxyPfpColl_t const &pfp_proxy = proxy::getCollection<std::vector<recob::PFParticle>>(evt, fPFPproducer,
+      proxy::withAssociated<larpandoraobj::PFParticleMetadata>(fPFPproducer),
+      proxy::withAssociated<recob::Cluster>(fCLSproducer),
+      proxy::withAssociated<recob::Slice>(fSLCproducer),
+      proxy::withAssociated<recob::Track>(fTRKproducer),
+      proxy::withAssociated<recob::Vertex>(fVTXproducer),
+      proxy::withAssociated<recob::PCAxis>(fPCAproducer),
+      proxy::withAssociated<recob::Shower>(fSHRproducer),
+      proxy::withAssociated<recob::SpacePoint>(fPFPproducer));
+    selection::BuildPFPMap(pfp_proxy);
+    // loopthrough PFParticles
+    for (const ProxyPfpElem_t &pfp_pxy : pfp_proxy)
+    {
+      // get metadata for this PFP
+      const auto &pfParticleMetadataList = pfp_pxy.get<larpandoraobj::PFParticleMetadata>();
+
+      //  find neutrino candidate
+      if (pfp_pxy->IsPrimary() == false) continue;
+      auto PDG = fabs(pfp_pxy->PdgCode());
+      fData->nu_pfp_pdg=PDG;
+      if ( (PDG == 12) || (PDG == 14) )
+      {
+        //std::cout<<"Found a neutrino PFP\n";
+        if (pfParticleMetadataList.size() != 0)
+        {
+          for (unsigned int j = 0; j < pfParticleMetadataList.size(); ++j)
+          {
+            const art::Ptr<larpandoraobj::PFParticleMetadata> &pfParticleMetadata(pfParticleMetadataList.at(j));
+            auto pfParticlePropertiesMap = pfParticleMetadata->GetPropertiesMap();
+            if (!pfParticlePropertiesMap.empty())
+            {
+              for (std::map<std::string, float>::const_iterator it = pfParticlePropertiesMap.begin(); it != pfParticlePropertiesMap.end(); ++it)
+              {
+                if( it->first == "IsNeutrino" ) fData->nu_isNeutrino  = it->second;
+                if( it->first == "NuScore"    ) fData->nu_nuscore.push_back(it->second);
+              }
+            }
+          }
+        } // if PFP metadata exists!
+  
+
+        // Get vertex info
+        double xyz[3] = {};
+        auto vtx = pfp_pxy.get<recob::Vertex>();
+        if (vtx.size() == 1)
+        {
+          // save vertex to array
+          vtx.at(0)->XYZ(xyz);
+          auto nuvtx = TVector3(xyz[0], xyz[1], xyz[2]);
+          float _reco_nu_vtx_sce[3];
+          nuselection::ApplySCECorrectionXYZ(nuvtx.X(),nuvtx.Y(),nuvtx.Z(), _reco_nu_vtx_sce);
+          fData->nu_reco_vtx_x = _reco_nu_vtx_sce[0];
+          fData->nu_reco_vtx_y = _reco_nu_vtx_sce[1];
+          fData->nu_reco_vtx_z = _reco_nu_vtx_sce[2];
+          //std::cout<<"Vertex:  "<<nuvtx.X()<<"  "<<nuvtx.Y()<<"  "<<nuvtx.Z()<<"\n";
+        }
+        else
+        {
+          std::cout << "ERROR. Found neutrino PFP w/ != 1 associated vertices..." << std::endl;
+        }
+
+        // collect PFParticle hierarchy originating from this neutrino candidate
+        std::vector<ProxyPfpElem_t> slice_pfp_v;
+        selection::AddDaughters(pfp_pxy, pfp_proxy, slice_pfp_v);
+        //std::cout << "This slice has " << slice_pfp_v.size() << " daughter PFParticles" << std::endl;
+        // create list of tracks and showers associated to this slice
+        std::vector<float>  trkscore_v;
+        std::vector<int>    trkid_v;
+        std::vector<float>  shwrscore_v;
+
+        for (auto pfp : slice_pfp_v)
+        { 
+          auto const &ass_trk_v = pfp.get<recob::Track>();
+          auto const &ass_shr_v = pfp.get<recob::Shower>();
+          float score = nuselection::GetTrackShowerScore(pfp);
+          if (ass_trk_v.size() == 1) {
+            trkid_v   .push_back(ass_trk_v.at(0)->ID());
+            trkscore_v.push_back(score);
+          }
+          if (ass_shr_v.size() == 1) {
+            shwrscore_v.push_back(score);
+          }
+        } // for all PFParticles in the slice
+        //std::cout<<"  - "<<trkscore_v.size()<<" tracks\n";
+        //std::cout<<"  - "<<shwrscore_v.size()<<" showers\n";
+        for(size_t itrk = 0; itrk < trkscore_v.size(); itrk++){
+          fData->nu_trk_id    .push_back(trkid_v[itrk]);
+          fData->nu_trk_score .push_back(trkscore_v[itrk]);
+        }
+        for(size_t ishwr = 0; ishwr < shwrscore_v.size(); ishwr++){
+          fData->nu_shwr_score .push_back(shwrscore_v[ishwr]);
+        }
+      
+      }//if PDG of neutrino
+    }//end loop over PFPs
+  }
+  
   
   //===========================================
   // Check if blip objects were saved to the event;
diff --git a/ubreco/BlipReco/BlipAna_module.cc b/ubreco/BlipReco/BlipAna_module.cc
index 5aa2f868..614d2901 100644
--- a/ubreco/BlipReco/BlipAna_module.cc
+++ b/ubreco/BlipReco/BlipAna_module.cc
@@ -292,7 +292,7 @@ class BlipAnaTreeDataStruct
   int   blip_proxtrkid[kMaxBlips];    // index of nearest trk
   bool  blip_touchtrk[kMaxBlips];     // is blip touching track?
   int   blip_touchtrkid[kMaxBlips];   // track ID of touched track
-  bool  blip_incylinder[kMaxBlips];   // is blip within a cylinder near a track
+  //bool  blip_incylinder[kMaxBlips];   // is blip within a cylinder near a track
   int   blip_clustid[kNplanes][kMaxBlips];     // cluster ID per plane
 
   // --- Reconstructed neutrino slice information (Pandora) --
@@ -483,7 +483,7 @@ class BlipAnaTreeDataStruct
     FillWith(blip_proxtrkid,  -9);
     FillWith(blip_touchtrk,   false);
     FillWith(blip_touchtrkid,  -9);
-    FillWith(blip_incylinder, false);
+    //FillWith(blip_incylinder, false);
     FillWith(blip_edepid,     -9);
     FillWith(blip_g4id,     -9);
     for(int i=0; i<kNplanes; i++){ 
@@ -609,7 +609,7 @@ class BlipAnaTreeDataStruct
     evtTree->Branch("blip_energy",blip_energy,"blip_energy[nblips]/F");
     evtTree->Branch("blip_yzcorr",blip_yzcorr,"blip_yzcorr[nblips]/F");
     //evtTree->Branch("blip_energyTrue",blip_energyTrue,"blip_energyTrue[nblips]/F");
-    evtTree->Branch("blip_incylinder",blip_incylinder,"blip_incylinder[nblips]/O");
+    //evtTree->Branch("blip_incylinder",blip_incylinder,"blip_incylinder[nblips]/O");
     evtTree->Branch("blip_proxtrkdist",blip_proxtrkdist,"blip_proxtrkdist[nblips]/F");
     evtTree->Branch("blip_touchtrk",blip_touchtrk,"blip_touchtrk[nblips]/O");
     if( saveTrkInfo ) {
@@ -812,6 +812,10 @@ class BlipAna : public art::EDAnalyzer
   //std::vector<std::unique_ptr<::analysis::AnalysisToolBase>> _analysisToolsVec;
 
   // --- Histograms ---
+  
+  TH1D*   h_ncapture_gammaE;
+  TH1D*   h_ncapture_gammaEsum;
+
   TH1D*   h_part_process;
  
   TH1D*   h_nhits[kNplanes];
@@ -971,7 +975,10 @@ class BlipAna : public art::EDAnalyzer
 
     // MC histograms related to truth
     art::TFileDirectory dir_truth = dir_diag.mkdir("Truth");
-    
+   
+    h_ncapture_gammaE    = dir_truth.make<TH1D>("ncapture_gammaE",   "Gammas from neutron capture;True gamma energy [MeV]",300,0,15);
+    h_ncapture_gammaEsum = dir_truth.make<TH1D>("ncapture_gammaEsum","Gammas from neutron capture;Summed gamma energy [MeV]",300,0,15);
+
     h_part_process    = dir_truth.make<TH1D>("part_process","MCParticle->Process()",5,0,5);
     auto xa = h_part_process->GetXaxis();
     xa->SetBinLabel(1,"primary");
@@ -1462,6 +1469,8 @@ void BlipAna::analyze(const art::Event& evt)
   }
 
 
+  std::map< int, std::vector<int> > map_neutron_gammas;
+  std::map< int, std::vector<float> > map_neutron_gammaE;
 
   //====================================
   // Save MCParticle information
@@ -1479,7 +1488,7 @@ void BlipAna::analyze(const art::Event& evt)
       map_g4trkid_index[pPart->TrackId()] = i;
       total_depEnergy       += pinfo[i].depEnergy;
       total_depElectrons    += pinfo[i].depElectrons;
-      
+
       // Save to TTree object
       if(i<kMaxG4){
         
@@ -1529,19 +1538,53 @@ void BlipAna::analyze(const art::Event& evt)
         fData->part_depElectrons[i]    = pinfo[i].depElectrons;
         fData->part_isPrimary[i]       = pinfo[i].isPrimary;
         // check containment
-        bool startInAV  = BlipUtils::IsPointInAV(pPart->Vx(),pPart->Vy(),pPart->Vz());
-        bool endInAV    = BlipUtils::IsPointInAV(pPart->EndPosition()[0],pPart->EndPosition()[1],pPart->EndPosition()[2]);
+        bool startInAV  = BlipUtils::IsPointInAV(pPart->Vx(),pPart->Vy(),pPart->Vz(),2);
+        bool endInAV    = BlipUtils::IsPointInAV(pPart->EndPosition()[0],pPart->EndPosition()[1],pPart->EndPosition()[2],2);
         if( startInAV && endInAV ) fData->part_isContained[i] = true;
- 
         if( fDebugMode && pPart->Process() != "muIoni" ) {
           PrintParticleInfo(i);
           printed++;
         }
+        
+        //-----------------------------------------------
+        // Check contained neutron captures
+        if( pPart->PdgCode() == 2112 && fData->part_isContained[i] && pPart->EndProcess() == "nCapture" ) {
+          map_neutron_gammas[pPart->TrackId()].clear();
+        }
+        // Check for neutron capture products
+        if( pPart->PdgCode() == 22 && pPart->Process() == "nCapture" ) {
+          if( map_neutron_gammas.count(pPart->Mother()) ) {
+            map_neutron_gammas[pPart->Mother()].push_back(pPart->TrackId());
+            map_neutron_gammaE[pPart->Mother()].push_back(pPart->P()*1e3);
+          }
+        }
+        //-----------------------------------------------
+      
+      
       }
+        
     } // endloop over G4 particles
-  
+    
+
+
     if( fDebugMode ) std::cout<<"True total energy deposited: "<<total_depEnergy<<" MeV \n";
-  
+    int nNeutrons = map_neutron_gammas.size();
+    //std::cout<<"Neutrons in this event: "<<map_neutron_gammas.size()<<"\n";
+    if( nNeutrons > 0 ) {
+      for(auto cascade : map_neutron_gammaE ) {
+        //std::cout<<" - "<<cascade.first<<": ";
+        float sumE = 0;
+        for(auto gammaE : cascade.second ) {
+          sumE += gammaE;
+          h_ncapture_gammaE->Fill(gammaE);
+          //std::cout<<gammaE<<" ";
+        }
+        h_ncapture_gammaEsum->Fill(sumE);
+        //std::cout<<"SumE = "<<sumE<<"\n";
+        //std::cout<<"\n";
+      }
+    }
+
   }//endif particles found in event
   
 
@@ -2162,7 +2205,7 @@ void BlipAna::analyze(const art::Event& evt)
     fData->blip_proxtrkid[i]  = blp.ProxTrkID;
     fData->blip_touchtrk[i]   = (blp.TouchTrkID >= 0 );
     fData->blip_touchtrkid[i] = blp.TouchTrkID;
-    fData->blip_incylinder[i] = blp.inCylinder;
+    //fData->blip_incylinder[i] = blp.inCylinder;
     fData->blip_charge[i]     = blp.Charge;
     fData->blip_energy[i]     = blp.Energy;
     fData->blip_yzcorr[i]     = tpcCalib.YZdqdxCorrection(fCaloPlane,blp.Position.Y(),blp.Position.Z());
diff --git a/ubreco/BlipReco/CMakeLists.txt b/ubreco/BlipReco/CMakeLists.txt
index 4cb531bb..d8e128cb 100644
--- a/ubreco/BlipReco/CMakeLists.txt
+++ b/ubreco/BlipReco/CMakeLists.txt
@@ -82,6 +82,7 @@ cet_build_plugin(
 
 add_subdirectory(Utils)
 add_subdirectory(Alg)
+add_subdirectory(RNN)
 add_subdirectory(job)
 install_fhicl()
 install_source()
diff --git a/ubreco/BlipReco/RNN/BlipRNNAlg.cc b/ubreco/BlipReco/RNN/BlipRNNAlg.cc
new file mode 100644
index 00000000..f3fd8d40
--- /dev/null
+++ b/ubreco/BlipReco/RNN/BlipRNNAlg.cc
@@ -0,0 +1,166 @@
+#include "ubreco/BlipReco/RNN/BlipRNNAlg.h"
+
+namespace blip {
+
+  //###########################################################
+  // Constructor
+  //###########################################################
+  BlipRNNAlg::BlipRNNAlg( fhicl::ParameterSet const& pset ):
+    device(torch::kCPU)
+  {
+    fModelFilename   = pset.get<std::string>   ("ModelFilename", "model_blipRNN.pt");
+    model_loaded=false;
+    loadModel();
+  }
+  
+  //====================================================
+  // Destructor
+  BlipRNNAlg::~BlipRNNAlg()
+  {
+  }
+  
+  //====================================================
+  // Load model
+  bool BlipRNNAlg::loadModel() {
+    // Skip if already loaded with same path
+    if (model_loaded && loaded_model_path == fModelFilename ) {
+      std::cout << "Model already loaded, skipping...\n";
+      return true;
+    }
+    try {
+      cet::search_path sp("FW_SEARCH_PATH");
+      std::cout<<"Searching for file...\n";
+      std::cout<<sp.find_file(fModelFilename)<<"\n";
+      model = torch::jit::load(sp.find_file(fModelFilename));
+      model.to(device);
+      model.eval();
+      model_loaded = true;
+      loaded_model_path = fModelFilename;
+      std::cout << "Model loaded: " << loaded_model_path << std::endl;
+      return true;
+    } catch (...) {
+        std::cout << "Failed to load model "<<fModelFilename<<"\n";
+        return false;
+      }
+  }
+  
+  //====================================================
+  // Predict blip 3D direction
+  std::vector<float> BlipRNNAlg::predict(vvfloat_t const& blip_hit_seq){
+    if( !model_loaded ) return {};
+    try{
+      
+      int nFeatures=5;
+      int nHits=blip_hit_seq.size();
+      auto bliphits = torch::zeros({1, nHits, nFeatures});
+      for(int i=0; i<nHits; i++){
+        for(int j=0; j<nFeatures; j++) 
+          bliphits[0][i][j] = blip_hit_seq[i][j];
+      }
+     
+      // create input
+      auto length = torch::tensor({nHits}, torch::kLong);
+      torch::NoGradGuard no_grad;
+      std::vector<torch::jit::IValue> input;
+      input.push_back(bliphits);
+      input.push_back(length);
+      
+      // feed it to the model
+      at::Tensor output = model.forward(input).toTensor();
+
+      // get the result and normalize magnitude
+      auto accessor = output.accessor<float, 2>();
+      TVector3 vdir(accessor[0][0],accessor[0][1],accessor[0][2]);
+      vdir.SetMag(1.);
+      
+      return { (float)vdir.X(), (float)vdir.Y(), (float)vdir.Z() };
+    
+    } catch (const std::exception& e) {
+      std::cerr << "RNN Prediction failed: "<<e.what()<<std::endl;
+      return {};
+    }
+  }
+
+  //============================================================
+  // This function takes in a blipobj::Blip object and pre-processed the
+  // hit information to feed into the function above
+  std::vector<float> BlipRNNAlg::predict(blipobj::Blip const& blp, std::vector<blipobj::HitInfo> const& hitinfo){
+    
+    // 3D direction (initialized to dummy)
+    std::vector<float> dir{-9,-9,-9};
+
+    // Basic cuts: require blip be valid
+    if( !blp.isValid ) return dir;
+
+    // Conversion factor for casting ticks into units of
+    // wire spacings (TODO: make this not hard-coded?)
+    float tick_to_us  = 0.5; // microsecond
+    float drift_speed = 1.1; // mm per microsecond
+    float wire_sep    = 3.0; // mm
+    float convfact    = tick_to_us * drift_speed / wire_sep;
+    
+    // Collections of hit info info per plane
+    std::array<std::vector<float>, kNplanes > blip_pl_wire;   // [plane] -> vector of wires
+    std::array<std::vector<float>, kNplanes > blip_pl_time;
+    std::array<std::vector<float>, kNplanes > blip_pl_rms;
+    std::array<std::vector<float>, kNplanes > blip_pl_amp;
+
+    // Loop over the planes
+    int nhits = 0;
+    float sum_driftT = 0;
+    float min_driftT = 999999;
+    for(size_t i=0; i<kNplanes; i++){
+      auto& clust = blp.clusters[i];
+      
+      // loop this blip's hits on this plane
+      for(auto hid : clust.HitIDs ){
+        float t = hitinfo[hid].driftTime * convfact;
+        float rms = hitinfo[hid].rms * convfact;
+        float amp = hitinfo[hid].amp;
+        blip_pl_wire[i].push_back( hitinfo[hid].wire );
+        blip_pl_time[i].push_back( t );
+        blip_pl_rms[i] .push_back( rms );
+        blip_pl_amp[i] .push_back( amp );
+        nhits++;
+        sum_driftT += t;
+        if( t < min_driftT ) min_driftT = t;
+      }
+      
+      // basic sanity check
+      if( blip_pl_wire[i].empty() ) continue;
+
+      // make corrections to wires on each plane
+      float max_wire = *std::max_element(blip_pl_wire[i].begin(), blip_pl_wire[i].end());
+      for(auto &wn : blip_pl_wire[i]) wn = max_wire - wn;
+    }
+
+    // If there were 2 hits or less, forget about it
+    if( nhits <= 2 ) return dir;
+    
+    // Ok we got all the hits. Shift the drift times
+    for(size_t i=0; i<kNplanes; i++){
+      for(auto& t : blip_pl_time[i] ) t -= min_driftT;
+    }
+    
+    // create vector of hit info needed for RNN inference
+    // order of hit features: rms, amp, wire, plane, time
+    vvfloat_t blip_hit_seq;
+    for(size_t i=0; i<kNplanes; i++){
+      for(size_t j=0; j<blip_pl_wire[i].size(); j++){
+        std::vector<float> hit_features;
+        hit_features.push_back( blip_pl_rms[i][j] );
+        hit_features.push_back( blip_pl_amp[i][j] );
+        hit_features.push_back( blip_pl_wire[i][j] );
+        hit_features.push_back( i );
+        hit_features.push_back( blip_pl_time[i][j] );
+        blip_hit_seq.push_back( hit_features );
+      }
+    }
+    
+    return predict(blip_hit_seq);
+  
+  }
+
+}
+
+
diff --git a/ubreco/BlipReco/RNN/BlipRNNAlg.h b/ubreco/BlipReco/RNN/BlipRNNAlg.h
new file mode 100644
index 00000000..a5f231e4
--- /dev/null
+++ b/ubreco/BlipReco/RNN/BlipRNNAlg.h
@@ -0,0 +1,48 @@
+#ifndef BLIPRNNALG_H
+#define BLIPRNNALG_H
+
+// Framework, libtorch, etc
+#include "fhiclcpp/ParameterSet.h"
+#include <torch/script.h>
+#include <algorithm>
+
+// Blip data types
+#include "ubreco/BlipReco/Utils/BlipCore.h"
+
+
+typedef std::vector<float>      vfloat_t;
+typedef std::vector<vfloat_t>   vvfloat_t;
+
+namespace blip {
+
+  //--------------------------------------------
+  class BlipRNNAlg {
+    
+    public:
+    
+      //Constructor/destructor
+      BlipRNNAlg( fhicl::ParameterSet const& pset );
+      BlipRNNAlg();
+      ~BlipRNNAlg();
+      
+      bool loadModel();
+      vfloat_t predict(vvfloat_t const&);
+      vfloat_t predict(blipobj::Blip const&, std::vector<blipobj::HitInfo> const&);
+
+      bool model_loaded;
+      std::string loaded_model_path;
+
+    private:
+     
+      // FCL configs
+      std::string fModelFilename;
+     
+      // PyTorch things
+      torch::jit::script::Module model;
+      torch::Device device;
+  
+  };
+
+}
+
+#endif
diff --git a/ubreco/BlipReco/RNN/CMakeLists.txt b/ubreco/BlipReco/RNN/CMakeLists.txt
new file mode 100644
index 00000000..093902bc
--- /dev/null
+++ b/ubreco/BlipReco/RNN/CMakeLists.txt
@@ -0,0 +1,17 @@
+cet_make_library(
+  SOURCE BlipRNNAlg.cc
+  LIBRARIES
+    PUBLIC
+    art_root_io::TFileService_service
+    ubobj::DataOverlay
+    lardata::Utilities
+    nusimdata::SimulationBase
+    TorchScatter::TorchScatter
+)
+
+install_headers()
+install_fhicl()
+install_source()
+
+file(GLOB model_file *.pt)
+install_fw(LIST ${model_file})
diff --git a/ubreco/BlipReco/RNN/model_blipRNN.pt b/ubreco/BlipReco/RNN/model_blipRNN.pt
new file mode 100644
index 00000000..01b1888d
Binary files /dev/null and b/ubreco/BlipReco/RNN/model_blipRNN.pt differ
diff --git a/ubreco/BlipReco/Utils/BlipCore.h b/ubreco/BlipReco/Utils/BlipCore.h
new file mode 100644
index 00000000..b6e029ec
--- /dev/null
+++ b/ubreco/BlipReco/Utils/BlipCore.h
@@ -0,0 +1,12 @@
+///////////////////////////////////////////////////////
+// BlipCore.h
+//
+// include-guard for 'DataTypes.h'
+//
+//////////////////////////////////////////////////////
+#ifndef BLIPCORE_H_SEEN
+#define BLIPCORE_H_SEEN
+
+#include "ubobj/Blip/DataTypes.h"
+
+#endif
diff --git a/ubreco/BlipReco/Utils/BlipUtils.cc b/ubreco/BlipReco/Utils/BlipUtils.cc
index a6e97a56..8965e61e 100644
--- a/ubreco/BlipReco/Utils/BlipUtils.cc
+++ b/ubreco/BlipReco/Utils/BlipUtils.cc
@@ -147,7 +147,6 @@ namespace BlipUtils {
       float w2 = pinfo.depEnergy/totE;
       tblip.Position    = w1*tblip.Position + w2*pinfo.position;
       tblip.Time        = w1*tblip.Time     + w2*pinfo.time;
-      tblip.LeadCharge  = pinfo.depElectrons;
     // ... if the particle isn't a match, show's over
     } else {
       return;
@@ -158,8 +157,7 @@ namespace BlipUtils {
     tblip.NumElectrons+= std::max(0.,pinfo.numElectrons);
 
     auto const detProp   = art::ServiceHandle<detinfo::DetectorPropertiesService>()->DataForJob();
-    float driftVel = detProp.DriftVelocity(detProp.Efield(0),detProp.Temperature());
-    tblip.DriftTime = tblip.Position.X() / driftVel;
+    tblip.DriftTime = tblip.Position.X() / detProp.DriftVelocity(detProp.Efield(0),detProp.Temperature());
 
     tblip.G4ChargeMap[part.TrackId()] += pinfo.depElectrons;
     if(pinfo.depElectrons > tblip.LeadCharge ) {
@@ -434,12 +432,42 @@ namespace BlipUtils {
     return newblip;
     
   }
+ 
+  MCP_t ReturnMCParticle(int particleID){
+    art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
+    return pi_serv->TrackIdToParticle(particleID);
+  }
+
+  MCPVec_t ReturnAllAncestors(int particleID) {
+    MCPVec_t out;
+    art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
+    const sim::ParticleList& plist = pi_serv->ParticleList();
+    if( !plist.HasParticle(particleID) )  return out;
+    auto p = pi_serv->TrackIdToParticle(particleID);
+    while( particleID > 0 ){
+      if( !plist.HasParticle(p.Mother() ) ) { return out; }
+      p = pi_serv->TrackIdToParticle(p.Mother());
+      out.push_back(p);
+    }
+    return out;
+  }
+ 
+  /*
+  MCPVec_t ReturnAllDaughters(int particleID) {
+    MCPVec_t out;
+    art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
+    const sim::ParticleList& plist = pi_serv->ParticleList();
+    // under construction
+    for(
+    return out; 
+  }
+  */
 
 
   //====================================================================
   // 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){
     art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
     const sim::ParticleList& plist = pi_serv->ParticleList();
     if( particleID == ancestorID  )       return true;
@@ -671,12 +699,19 @@ namespace BlipUtils {
   }
 
   //===========================================================================
-  bool IsPointInAV(float x, float y, float z){
+  bool IsPointInAV(float x, float y, float z, float margin){
     
     // Get geo boundaries
     double xmin, xmax, ymin, ymax, zmin, zmax;
     GetGeoBoundaries(xmin,xmax,ymin,ymax,zmin,zmax);
-      
+     
+    // add 2cm margin
+    if(margin){
+      xmin += margin; xmax -= margin;
+      ymin += margin; ymax -= margin;
+      zmin += margin; zmax -= margin;
+    }
+
     if(     x >= xmin && x <= xmax
         &&  y >= ymin && y <= ymax
         &&  z >= zmin && z <= zmax ) {
@@ -687,8 +722,8 @@ namespace BlipUtils {
     
   }
   
-  bool IsPointInAV(TVector3& v){
-    return IsPointInAV(v.X(), v.Y(), v.Z());
+  bool IsPointInAV(TVector3& v,float margin){
+    return IsPointInAV(v.X(), v.Y(), v.Z(), margin);
   }
   
   
diff --git a/ubreco/BlipReco/Utils/BlipUtils.h b/ubreco/BlipReco/Utils/BlipUtils.h
index 612cfc04..5e7afba5 100644
--- a/ubreco/BlipReco/Utils/BlipUtils.h
+++ b/ubreco/BlipReco/Utils/BlipUtils.h
@@ -30,11 +30,16 @@
 #include <vector>
 #include <map>
 
-#include "ubobj/Blip/DataTypes.h"
+// BlipReco data types
+#include "ubreco/BlipReco/Utils/BlipCore.h"
+
+// ROOT
 #include "TH1D.h"
 
 
 typedef std::vector<art::Ptr<sim::SimEnergyDeposit>> SEDVec_t;
+typedef std::vector<simb::MCParticle> MCPVec_t;
+typedef simb::MCParticle MCP_t;
 
 namespace BlipUtils {
  
@@ -68,19 +73,22 @@ 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 breakAtPhots=false);
   double  DistToBoundary(const recob::Track::Point_t&);
   double  DistToLine(TVector3&, TVector3&, TVector3&);
   double  DistToLine2D(TVector2&, TVector2&, TVector2&);
   void    GetGeoBoundaries(double&,double&,double&,double&,double&,double&);
-  bool    IsPointInAV(float,float,float);
-  bool    IsPointInAV(TVector3&);
+  bool    IsPointInAV(float x,float y,float z,float margin=0);
+  bool    IsPointInAV(TVector3& p,float margin=0);
   bool    IsPointAtBnd(float,float,float);
   bool    IsPointAtBnd(TVector3&);
   void    NormalizeHist(TH1D*);
   float   FindMedian(std::vector<float>&);
   float   FindMean(std::vector<float>&);
-  
+  MCP_t   ReturnMCParticle(int);
+
+  MCPVec_t ReturnAllAncestors(int);
+  //MCPVec_t ReturnAllDaughters(int); 
 }
 
 #endif
diff --git a/ubreco/BlipReco/Utils/NuSelectionToolBase.h b/ubreco/BlipReco/Utils/NuSelectionToolBase.h
index 724ccc6b..b34bb6b0 100755
--- a/ubreco/BlipReco/Utils/NuSelectionToolBase.h
+++ b/ubreco/BlipReco/Utils/NuSelectionToolBase.h
@@ -33,6 +33,15 @@ namespace selection {
   using ProxyPfpElem_t = nuselection::ProxyPfpElem_t;
   using ProxyClusColl_t = nuselection::ProxyClusColl_t;
   using ProxyClusElem_t = nuselection::ProxyClusElem_t;
+  
+  // a map linking the PFP Self() attribute used for hierarchy building to the PFP index in the event record
+  std::map<unsigned int, unsigned int> _pfpmap;
+  
+  void BuildPFPMap(const ProxyPfpColl_t &pfp_pxy_col);
+  void AddDaughters(const ProxyPfpElem_t &pfp_pxy,
+                    const ProxyPfpColl_t &pfp_pxy_col,
+                    std::vector<ProxyPfpElem_t> &slice_v);
+
 
 class SelectionToolBase {
 
@@ -84,6 +93,42 @@ class SelectionToolBase {
 
 };
 
+
+
+void BuildPFPMap(const ProxyPfpColl_t &pfp_pxy_col)
+{
+  _pfpmap.clear();
+  unsigned int p = 0;
+  for (const auto &pfp_pxy : pfp_pxy_col)
+  { _pfpmap[pfp_pxy->Self()] = p; p++; }
+  return;
+} // BuildPFPMap
+
+void AddDaughters(const ProxyPfpElem_t &pfp_pxy,
+                  const ProxyPfpColl_t &pfp_pxy_col,
+                  std::vector<ProxyPfpElem_t> &slice_v)
+{
+  auto daughters = pfp_pxy->Daughters();
+  slice_v.push_back(pfp_pxy);
+  //std::cout << "\t PFP w/ PdgCode " << pfp_pxy->PdgCode() << " has " << daughters.size() << " daughters" << std::endl;
+  for (auto const &daughterid : daughters){
+    if (_pfpmap.find(daughterid) == _pfpmap.end())
+    {
+      //std::cout << "Did not find DAUGHTERID in map! error" << std::endl;
+      continue;
+    }
+
+    // const art::Ptr<recob::PFParticle> pfp_pxy(pfp_pxy_col, _pfpmap.at(daughterid) );
+    auto pfp_pxy2 = pfp_pxy_col.begin();
+    for (size_t j = 0; j < _pfpmap.at(daughterid); ++j)
+      ++pfp_pxy2;
+    // const T& pfp_pxy2 = (pfp_pxy_col.begin()+_pfpmap.at(daughterid));
+    AddDaughters(*pfp_pxy2, pfp_pxy_col, slice_v);
+  } // for all daughters
+
+  return;
+} // AddDaughters
+
 } // selection namespace
 
 #endif
diff --git a/ubreco/BlipReco/microboone_blipreco.fcl b/ubreco/BlipReco/microboone_blipreco.fcl
index 5565cd5b..f434b862 100644
--- a/ubreco/BlipReco/microboone_blipreco.fcl
+++ b/ubreco/BlipReco/microboone_blipreco.fcl
@@ -52,8 +52,8 @@ microboone_blipalg:
   ClustMatchQDiffLimit: 10e3            #// if charge-diff >, clusters are subject to QRatio cut above [e-]
   MinMatchedPlanes:     2               #// minimum number of planes matched for 3D blip
   PickyBlips:           false           #// require 3 planes that all intersect at ~same point
-  ApplyTrkCylinderCut:  false           #// reject blips in track cone-cylinder regions
-  CylinderRadius:       15.0            #// 3D cylinder cut radius [cm]
+  #ApplyTrkCylinderCut:  false           #// reject blips in track cone-cylinder regions
+  #CylinderRadius:       15.0            #// 3D cylinder cut radius [cm]
   IgnoreDataTrks:       false           #// only consider MC tracks when determining blip-track distances
   CaloPlane:            2               #// plane used for calorimetry (2=collection)
   CalodEdx:             2.8             #// assumed dE/dx for MeV-scale recombination corrections