diff --git a/CMakeLists.txt b/CMakeLists.txt
index 4c86b5b80..fcf0f5774 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -59,6 +59,7 @@ find_package( lardata REQUIRED )
 find_package( larevt REQUIRED )
 find_package( larsim REQUIRED )
 find_package( larreco REQUIRED )
+find_package( larrecodnn REQUIRED )
 find_package( larpandora REQUIRED )
 find_package( larpandoracontent REQUIRED )
 find_package( artdaq_core EXPORT REQUIRED )
diff --git a/sbncode/CAFMaker/CAFMakerParams.h b/sbncode/CAFMaker/CAFMakerParams.h
index f5d049d97..271e54353 100644
--- a/sbncode/CAFMaker/CAFMakerParams.h
+++ b/sbncode/CAFMaker/CAFMakerParams.h
@@ -134,6 +134,12 @@ namespace caf
       "pandora"
     };
 
+    Atom<string> CNNScoreLabel {
+      Name("CNNScoreLabel"),
+      Comment("Base label of CNN score producer."),
+      "cnnid"
+    };
+
     Atom<string> StubLabel {
       Name("StubLabel"),
       Comment("Base label of Stub producer."),
@@ -252,6 +258,18 @@ namespace caf
       "pandoraTrackCRTTrack"
     };
 
+    Atom<string> CRTSpacePointMatchLabel {
+      Name("CRTSpacePointMatchLabel"),
+      Comment("Base label of track to CRT spacepoint matching producer."),
+      "crtspacepointmatching"
+    };
+
+    Atom<string> SBNDCRTTrackMatchLabel {
+      Name("SBNDCRTTrackMatchLabel"),
+      Comment("Base label of track to SBND CRT track matching producer."),
+      "crttrackmatching"
+    };
+
     Atom<string> TrackMCSLabel {
       Name("TrackMCSLabel"),
       Comment("Base label of track MCS momentum calculation producer."),
@@ -267,15 +285,27 @@ namespace caf
     Atom<string> CRTHitLabel {
       Name("CRTHitLabel"),
       Comment("Label of sbn CRT hits."),
-      "crthit" // same for icarus and sbnd
+      "crthit" // icarus
     };
 
     Atom<string> CRTTrackLabel {
       Name("CRTTrackLabel"),
       Comment("Label of sbn CRT tracks."),
-      "crttrack" // same for icarus and sbnd
+      "crttrack" // icarus
     };
-    
+
+    Atom<string> CRTSpacePointLabel {
+      Name("CRTSpacePointLabel"),
+      Comment("Label of sbnd CRT spacepoints."),
+      "crtspacepoints" // sbnd
+    };
+
+    Atom<string> SBNDCRTTrackLabel {
+      Name("SBNDCRTTrackLabel"),
+      Comment("Label of sbnd CRT tracks."),
+      "crttracks" // sbnd
+    };
+
     Atom<string> CRTPMTLabel {
       Name("CRTPMTLabel"),
       Comment("Label for the CRTPMT Matched variables from the crtpmt data product"),
diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc
index c902272fc..a5b23527d 100644
--- a/sbncode/CAFMaker/CAFMaker_module.cc
+++ b/sbncode/CAFMaker/CAFMaker_module.cc
@@ -106,6 +106,7 @@
 #include "sbnobj/Common/SBNEventWeight/EventWeightMap.h"
 #include "sbnobj/Common/SBNEventWeight/EventWeightParameterSet.h"
 #include "sbnobj/Common/Reco/MVAPID.h"
+#include "sbnobj/Common/Reco/CNNScore.h"
 #include "sbnobj/Common/Reco/ScatterClosestApproach.h"
 #include "sbnobj/Common/Reco/StoppingChi2Fit.h"
 #include "sbnobj/Common/POTAccounting/BNBSpillInfo.h"
@@ -276,6 +277,10 @@ class CAFMaker : public art::EDProducer {
   art::FindOneP<T> FindOnePStrict(const U& from, const art::Event& evt,
 				  const art::InputTag& label) const;
 
+  template <class T, class D, class U>
+  art::FindOneP<T, D> FindOnePDStrict(const U& from,
+                                      const art::Event& evt,
+                                      const art::InputTag& tag) const;
 
   /// \brief Retrieve an object from an association, with error handling
   ///
@@ -1005,6 +1010,25 @@ art::FindOneP<T> CAFMaker::FindOnePStrict(const U& from,
   return ret;
 }
 
+//......................................................................
+template <class T, class D, class U>
+art::FindOneP<T, D> CAFMaker::FindOnePDStrict(const U& from,
+                                              const art::Event& evt,
+                                              const art::InputTag& tag) const {
+  art::FindOneP<T, D> ret(from, evt, tag);
+
+  if (!tag.label().empty() && !ret.isValid() && fParams.StrictMode()) {
+    std::cout << "CAFMaker: No Assn from '"
+              << cet::demangle_symbol(typeid(from).name()) << "' to '"
+              << cet::demangle_symbol(typeid(T).name())
+              << "' found under label '" << tag << "'. "
+              << "Set 'StrictMode: false' to continue anyway." << std::endl;
+    abort();
+  }
+
+  return ret;
+}
+
 //......................................................................
 template <class T>
 bool CAFMaker::GetAssociatedProduct(const art::FindManyP<T>& fm, int idx,
@@ -1311,40 +1335,64 @@ void CAFMaker::produce(art::Event& evt) noexcept {
     pass_flash_trig = *flashtrig_handle;
   }
 
-  // Fill various detector information associated with the event
-  //
-  // Get all of the CRT hits
-  std::vector<caf::SRCRTHit> srcrthits;
-
-  art::Handle<std::vector<sbn::crt::CRTHit>> crthits_handle;
-  GetByLabelStrict(evt, fParams.CRTHitLabel(), crthits_handle);
-  // fill into event
-  //int64_t CRT_T0_reference_time = fParams.ReferenceCRTT0ToBeam() ? -srtrigger.beam_gate_time_abs : 0; // ns, signed
-  //double CRT_T1_reference_time = fParams.ReferenceCRTT1FromTriggerToBeam() ? srtrigger.trigger_within_gate : 0.;
   int64_t CRT_T0_reference_time = isRealData ?  -srtrigger.beam_gate_time_abs : -fParams.CRTSimT0Offset();
   double CRT_T1_reference_time = isRealData ? srtrigger.trigger_within_gate : -fParams.CRTSimT0Offset();
 
-  if (crthits_handle.isValid()) {
-    const std::vector<sbn::crt::CRTHit> &crthits = *crthits_handle;
-    for (unsigned i = 0; i < crthits.size(); i++) {
-      srcrthits.emplace_back();
-      FillCRTHit(crthits[i], fParams.CRTUseTS0(), CRT_T0_reference_time, CRT_T1_reference_time, srcrthits.back());
-    }
-  }
+  // Fill various detector information associated with the event
 
-  // Get all of the CRT Tracks
+  std::vector<caf::SRCRTHit> srcrthits;
   std::vector<caf::SRCRTTrack> srcrttracks;
+  std::vector<caf::SRCRTSpacePoint> srcrtspacepoints;
+  std::vector<caf::SRSBNDCRTTrack> srsbndcrttracks;
 
-  art::Handle<std::vector<sbn::crt::CRTTrack>> crttracks_handle;
-  GetByLabelStrict(evt, fParams.CRTTrackLabel(), crttracks_handle);
-  // fill into event
-  if (crttracks_handle.isValid()) {
-    const std::vector<sbn::crt::CRTTrack> &crttracks = *crttracks_handle;
-    for (unsigned i = 0; i < crttracks.size(); i++) {
-      srcrttracks.emplace_back();
-      FillCRTTrack(crttracks[i], fParams.CRTUseTS0(), srcrttracks.back());
+  if(fDet == kICARUS)
+    {
+      art::Handle<std::vector<sbn::crt::CRTHit>> crthits_handle;
+      GetByLabelStrict(evt, fParams.CRTHitLabel(), crthits_handle);
+      // fill into event
+      if (crthits_handle.isValid()) {
+        const std::vector<sbn::crt::CRTHit> &crthits = *crthits_handle;
+        for (unsigned i = 0; i < crthits.size(); i++) {
+          srcrthits.emplace_back();
+          FillCRTHit(crthits[i], fParams.CRTUseTS0(), CRT_T0_reference_time, CRT_T1_reference_time, srcrthits.back());
+        }
+      }
+
+      art::Handle<std::vector<sbn::crt::CRTTrack>> crttracks_handle;
+      GetByLabelStrict(evt, fParams.CRTTrackLabel(), crttracks_handle);
+      // fill into event
+      if (crttracks_handle.isValid()) {
+        const std::vector<sbn::crt::CRTTrack> &crttracks = *crttracks_handle;
+        for (unsigned i = 0; i < crttracks.size(); i++) {
+          srcrttracks.emplace_back();
+          FillCRTTrack(crttracks[i], fParams.CRTUseTS0(), srcrttracks.back());
+        }
+      }
+    }
+  else if(fDet == kSBND)
+    {
+      art::Handle<std::vector<sbnd::crt::CRTSpacePoint>> crtspacepoints_handle;
+      GetByLabelStrict(evt, fParams.CRTSpacePointLabel(), crtspacepoints_handle);
+
+      if (crtspacepoints_handle.isValid()) {
+        const std::vector<sbnd::crt::CRTSpacePoint> &crtspacepoints = *crtspacepoints_handle;
+        for (unsigned i = 0; i < crtspacepoints.size(); i++) {
+          srcrtspacepoints.emplace_back();
+          FillCRTSpacePoint(crtspacepoints[i], srcrtspacepoints.back());
+        }
+      }
+
+      art::Handle<std::vector<sbnd::crt::CRTTrack>> sbndcrttracks_handle;
+      GetByLabelStrict(evt, fParams.SBNDCRTTrackLabel(), sbndcrttracks_handle);
+      // fill into event
+      if (sbndcrttracks_handle.isValid()) {
+        const std::vector<sbnd::crt::CRTTrack> &sbndcrttracks = *sbndcrttracks_handle;
+        for (unsigned i = 0; i < sbndcrttracks.size(); i++) {
+          srsbndcrttracks.emplace_back();
+          FillSBNDCRTTrack(sbndcrttracks[i], srsbndcrttracks.back());
+        }
+      }
     }
-  }
 
   // Get all of the CRTPMT Matches .. 
   std::vector<caf::SRCRTPMTMatch> srcrtpmtmatches;
@@ -1599,6 +1647,10 @@ void CAFMaker::produce(art::Event& evt) noexcept {
       FindManyPStrict<sbn::MVAPID>(slcShowers, evt,
           fParams.ShowerRazzleLabel() + slice_tag_suff);
 
+    art::FindOneP<sbn::PFPCNNScore> foCNNScores = 
+      FindOnePStrict<sbn::PFPCNNScore>(fmPFPart, evt,
+          fParams.CNNScoreLabel() + slice_tag_suff);
+
     art::FindManyP<recob::Vertex> fmVertex =
       FindManyPStrict<recob::Vertex>(fmPFPart, evt,
              fParams.PFParticleLabel() + slice_tag_suff);
@@ -1622,6 +1674,14 @@ void CAFMaker::produce(art::Event& evt) noexcept {
       FindManyPStrict<anab::T0>(slcTracks, evt,
                fParams.CRTTrackMatchLabel() + slice_tag_suff);
 
+    art::FindOneP<sbnd::crt::CRTSpacePoint, anab::T0> foCRTSpacePointMatch =
+      FindOnePDStrict<sbnd::crt::CRTSpacePoint, anab::T0>(slcTracks, evt,
+               fParams.CRTSpacePointMatchLabel() + slice_tag_suff);
+
+    art::FindOneP<sbnd::crt::CRTTrack, anab::T0> foSBNDCRTTrackMatch =
+      FindOnePDStrict<sbnd::crt::CRTTrack, anab::T0>(slcTracks, evt,
+               fParams.SBNDCRTTrackMatchLabel() + slice_tag_suff);
+
     std::vector<art::FindManyP<recob::MCSFitResult>> fmMCSs;
     static const std::vector<std::string> PIDnames {"muon", "pion", "kaon", "proton"};
     for (std::string pid: PIDnames) {
@@ -1788,6 +1848,11 @@ void CAFMaker::produce(art::Event& evt) noexcept {
       const larpandoraobj::PFParticleMetadata *pfpMeta = (fmPFPMeta.at(iPart).empty()) ? NULL : fmPFPMeta.at(iPart).at(0).get();
       FillPFPVars(thisParticle, primary, pfpMeta, thisPFPT0, pfp);
 
+      if (foCNNScores.isValid()) {
+        const sbn::PFPCNNScore *cnnScores = foCNNScores.at(iPart).get();
+        FillCNNScores(thisParticle, cnnScores, pfp);
+      }
+    
       if (!thisTrack.empty())  { // it has a track!
         assert(thisTrack.size() == 1);
 
@@ -1836,7 +1901,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
               fParams.TrackHitFillRRStartCut(), fParams.TrackHitFillRREndCut(),
               lar::providerFrom<geo::Geometry>(), dprop, trk);
         }
-        if (fmCRTHitMatch.isValid()) {
+        if (fmCRTHitMatch.isValid() && fDet == kICARUS) {
           art::FindManyP<sbn::crt::CRTHit> CRTT02Hit = FindManyPStrict<sbn::crt::CRTHit>
               (fmCRTHitMatch.at(iPart), evt, fParams.CRTHitMatchLabel() + slice_tag_suff);
          
@@ -1846,9 +1911,25 @@ void CAFMaker::produce(art::Event& evt) noexcept {
           FillTrackCRTHit(fmCRTHitMatch.at(iPart), crthitmatch, fParams.CRTUseTS0(), CRT_T0_reference_time, CRT_T1_reference_time, trk);
         }
         // NOTE: SEE TODO AT fmCRTTrackMatch
-        if (fmCRTTrackMatch.isValid()) {
+        if (fmCRTTrackMatch.isValid() && fDet == kICARUS) {
           FillTrackCRTTrack(fmCRTTrackMatch.at(iPart), trk);
         }
+
+        if(foCRTSpacePointMatch.isValid() && fDet == kSBND)
+          {
+            const art::Ptr<sbnd::crt::CRTSpacePoint> crtspacepoint = foCRTSpacePointMatch.at(iPart);
+
+            if(crtspacepoint.isNonnull())
+              FillTrackCRTSpacePoint(foCRTSpacePointMatch.data(iPart).ref(), crtspacepoint, trk);
+          }
+        if(foSBNDCRTTrackMatch.isValid() && fDet == kSBND)
+          {
+            const art::Ptr<sbnd::crt::CRTTrack> sbndcrttrack = foSBNDCRTTrackMatch.at(iPart);
+
+            if(sbndcrttrack.isNonnull())
+              FillTrackSBNDCRTTrack(foSBNDCRTTrackMatch.data(iPart).ref(), sbndcrttrack, trk);
+          }
+
         // Truth matching
         if (fmTrackHit.isValid()) {
           if ( !isRealData ) {
@@ -1912,17 +1993,21 @@ void CAFMaker::produce(art::Event& evt) noexcept {
   //#######################################################
   //  Fill rec Tree
   //#######################################################
-  rec.nslc            = rec.slc.size();
-  rec.mc              = srtruthbranch;
-  rec.fake_reco       = srfakereco;
-  rec.nfake_reco      = srfakereco.size();
-  rec.pass_flashtrig  = pass_flash_trig;  // trigger result
-  rec.crt_hits        = srcrthits;
-  rec.ncrt_hits       = srcrthits.size();
-  rec.crt_tracks        = srcrttracks;
-  rec.ncrt_tracks       = srcrttracks.size();
-  rec.opflashes       = srflashes;
-  rec.nopflashes      = srflashes.size();
+  rec.nslc             = rec.slc.size();
+  rec.mc               = srtruthbranch;
+  rec.fake_reco        = srfakereco;
+  rec.nfake_reco       = srfakereco.size();
+  rec.pass_flashtrig   = pass_flash_trig;  // trigger result
+  rec.crt_hits         = srcrthits;
+  rec.ncrt_hits        = srcrthits.size();
+  rec.crt_tracks       = srcrttracks;
+  rec.ncrt_tracks      = srcrttracks.size();
+  rec.crt_spacepoints  = srcrtspacepoints;
+  rec.ncrt_spacepoints = srcrtspacepoints.size();
+  rec.sbnd_crt_tracks  = srsbndcrttracks;
+  rec.nsbnd_crt_tracks = srsbndcrttracks.size();
+  rec.opflashes        = srflashes;
+  rec.nopflashes       = srflashes.size();
   if (fParams.FillTrueParticles()) {
     rec.true_particles  = true_particles;
   }
diff --git a/sbncode/CAFMaker/CMakeLists.txt b/sbncode/CAFMaker/CMakeLists.txt
index cd0283055..de8932790 100644
--- a/sbncode/CAFMaker/CMakeLists.txt
+++ b/sbncode/CAFMaker/CMakeLists.txt
@@ -34,6 +34,7 @@ art_make_library( LIBRARY_NAME sbncode_CAFMaker
                   sbnobj::Common_CRT
                   sbnobj::Common_Reco
                   sbnobj::Common_Analysis
+                  sbnobj::SBND_CRT
                   lardataalg::DetectorInfo
                   art::Framework_Services_System_TriggerNamesService_service
                   sbncode_Metadata_MetadataSBN_service
diff --git a/sbncode/CAFMaker/FillReco.cxx b/sbncode/CAFMaker/FillReco.cxx
index 0655b1c02..6c085f655 100644
--- a/sbncode/CAFMaker/FillReco.cxx
+++ b/sbncode/CAFMaker/FillReco.cxx
@@ -112,6 +112,31 @@ namespace caf
     srtrack.hitb.plane = track.plane2;
   }
 
+  void FillCRTSpacePoint(const sbnd::crt::CRTSpacePoint &spacepoint,
+                         caf::SRCRTSpacePoint &srspacepoint,
+                         bool allowEmpty)
+  {
+    srspacepoint.position     = SRVector3D(spacepoint.X(), spacepoint.Y(), spacepoint.Z());
+    srspacepoint.position_err = SRVector3D(spacepoint.XErr(), spacepoint.YErr(), spacepoint.ZErr());
+    srspacepoint.pe           = spacepoint.PE();
+    srspacepoint.time         = spacepoint.Time();
+    srspacepoint.time_err     = spacepoint.TimeErr();
+    srspacepoint.complete     = spacepoint.Complete();
+  }
+
+  void FillSBNDCRTTrack(const sbnd::crt::CRTTrack &track,
+                        caf::SRSBNDCRTTrack &srsbndcrttrack,
+                        bool allowEmpty)
+  {
+    for(auto const& point : track.Points())
+      srsbndcrttrack.points.emplace_back(point.X(), point.Y(), point.Z());
+
+    srsbndcrttrack.time    = track.Time();
+    srsbndcrttrack.time_err = track.TimeErr();
+    srsbndcrttrack.pe       = track.PE();
+    srsbndcrttrack.tof      = track.ToF();
+  }
+
   void FillCRTPMTMatch(const sbn::crt::CRTPMTMatching &match,
 		       caf::SRCRTPMTMatch &srmatch,
 		       bool allowEmpty){
@@ -540,6 +565,24 @@ namespace caf
     }
   }
 
+  void FillTrackCRTSpacePoint(const anab::T0 &t0match,
+                              const art::Ptr<sbnd::crt::CRTSpacePoint> &spacepointmatch,
+                              caf::SRTrack &srtrack,
+                              bool allowEmpty)
+  {
+    srtrack.crtspacepoint.score = t0match.fTriggerConfidence;
+    FillCRTSpacePoint(*spacepointmatch, srtrack.crtspacepoint.spacepoint);
+  }
+
+  void FillTrackSBNDCRTTrack(const anab::T0 &t0match,
+                             const art::Ptr<sbnd::crt::CRTTrack> &trackmatch,
+                             caf::SRTrack &srtrack,
+                             bool allowEmpty)
+  {
+    srtrack.crtsbndtrack.score = t0match.fTriggerConfidence;
+    FillSBNDCRTTrack(*trackmatch, srtrack.crtsbndtrack.track);
+  }
+
   void FillTrackMCS(const recob::Track& track,
                     const std::array<std::vector<art::Ptr<recob::MCSFitResult>>, 4> &mcs_results,
                     caf::SRTrack& srtrack,
@@ -875,6 +918,24 @@ namespace caf
     }
   }
 
+  void FillCNNScores(const recob::PFParticle &particle,
+                     const sbn::PFPCNNScore *cnnscore,
+                     caf::SRPFP& srpfp,
+                     bool allowEmpty)
+  {
+    // std::cout << "pfpTrackScore: " << cnnscore->pfpTrackScore << std::endl;
+    // std::cout << "pfpShowerScore: " << cnnscore->pfpShowerScore << std::endl;
+    // std::cout << "pfpNoiseScore: " << cnnscore->pfpNoiseScore << std::endl;
+    // std::cout << "pfpMichelScore: " << cnnscore->pfpMichelScore << std::endl;
+
+    srpfp.cnnscore.track = cnnscore->pfpTrackScore;
+    srpfp.cnnscore.shower = cnnscore->pfpShowerScore;
+    srpfp.cnnscore.noise = cnnscore->pfpNoiseScore;
+    srpfp.cnnscore.michel = cnnscore->pfpMichelScore;
+    srpfp.cnnscore.endmichel = cnnscore->pfpEndMichelScore;
+    srpfp.cnnscore.nclusters = cnnscore->nClusters;
+  }
+
   void FillHitVars(const recob::Hit& hit,
                    unsigned producer,
                    const recob::SpacePoint& spacepoint,
diff --git a/sbncode/CAFMaker/FillReco.h b/sbncode/CAFMaker/FillReco.h
index ff68916a1..a90de4ada 100644
--- a/sbncode/CAFMaker/FillReco.h
+++ b/sbncode/CAFMaker/FillReco.h
@@ -30,6 +30,7 @@
 #include "sbnobj/Common/Reco/RangeP.h"
 #include "sbnobj/Common/Reco/ShowerSelectionVars.h"
 #include "sbnobj/Common/Reco/MVAPID.h"
+#include "sbnobj/Common/Reco/CNNScore.h"
 #include "sbnobj/Common/Reco/ScatterClosestApproach.h"
 #include "sbnobj/Common/Reco/StoppingChi2Fit.h"
 #include "sbnobj/Common/Reco/CRUMBSResult.h"
@@ -37,6 +38,8 @@
 #include "sbnobj/Common/Reco/TPCPMTBarycenterMatch.h"
 #include "sbnobj/Common/CRT/CRTHit.hh"
 #include "sbnobj/Common/CRT/CRTTrack.hh"
+#include "sbnobj/SBND/CRT/CRTSpacePoint.hh"
+#include "sbnobj/SBND/CRT/CRTTrack.hh"
 #include "sbnobj/Common/CRT/CRTPMTMatching.hh"
 #include "nusimdata/SimulationBase/MCParticle.h"
 #include "nusimdata/SimulationBase/MCTruth.h"
@@ -123,6 +126,11 @@ namespace caf
                    caf::SRPFP& srpfp,
                    bool allowEmpty= false);
 
+  void FillCNNScores(const recob::PFParticle &particle,
+                     const sbn::PFPCNNScore *cnnscore,
+                     caf::SRPFP& srpfp,
+                     bool allowEmpty = false);
+
   void FillTrackCRTHit(const std::vector<art::Ptr<anab::T0>> &t0match,
                        const std::vector<art::Ptr<sbn::crt::CRTHit>> &hitmatch,
                        bool use_ts0,
@@ -135,6 +143,16 @@ namespace caf
                        caf::SRTrack &srtrack,
                        bool allowEmpty = false);
 
+  void FillTrackCRTSpacePoint(const anab::T0 &t0match,
+                              const art::Ptr<sbnd::crt::CRTSpacePoint> &spacepointmatch,
+                              caf::SRTrack &srtrack,
+                              bool allowEmpty = false);
+
+  void FillTrackSBNDCRTTrack(const anab::T0 &t0match,
+                             const art::Ptr<sbnd::crt::CRTTrack> &trackmatch,
+                             caf::SRTrack &srtrack,
+                             bool allowEmpty = false);
+
   void FillTrackMCS(const recob::Track& track,
                     const std::array<std::vector<art::Ptr<recob::MCSFitResult>>, 4> &mcs_results,
                     caf::SRTrack& srtrack,
@@ -193,6 +211,14 @@ namespace caf
                   caf::SRCRTTrack &srtrack,
                   bool allowEmpty = false);
 
+  void FillCRTSpacePoint(const sbnd::crt::CRTSpacePoint &spacepoint,
+                         caf::SRCRTSpacePoint &srspacepoint,
+                         bool allowEmpty = false);
+
+  void FillSBNDCRTTrack(const sbnd::crt::CRTTrack &track,
+                        caf::SRSBNDCRTTrack &srsbndcrttrack,
+                        bool allowEmpty = false);
+
   void FillOpFlash(const recob::OpFlash &flash,
                   std::vector<recob::OpHit const*> const& hits,
                   int cryo,
diff --git a/sbncode/CosmicID/CMakeLists.txt b/sbncode/CosmicID/CMakeLists.txt
index 4ab97c17f..527721870 100644
--- a/sbncode/CosmicID/CMakeLists.txt
+++ b/sbncode/CosmicID/CMakeLists.txt
@@ -4,6 +4,7 @@ cet_build_plugin( CRUMBS art::module
           fhiclcpp::fhiclcpp
           art::Persistency_Provenance canvas::canvas
           sbnobj::Common_Reco
+          sbnobj::SBND_CRT
           sbncode_GeoWrappers
           sbncode_LArRecoProducer
           lardataobj::RecoBase_AssnsDicts_dict
diff --git a/sbncode/CosmicID/CRUMBS_module.cc b/sbncode/CosmicID/CRUMBS_module.cc
index 469787d6b..d4c583b34 100644
--- a/sbncode/CosmicID/CRUMBS_module.cc
+++ b/sbncode/CosmicID/CRUMBS_module.cc
@@ -40,6 +40,8 @@
 #include "sbnobj/Common/Reco/SimpleFlashMatchVars.h"
 #include "sbnobj/Common/Reco/StoppingChi2Fit.h"
 #include "sbncode/LArRecoProducer/TrackStoppingChi2Alg.h"
+#include "sbnobj/SBND/CRT/CRTTrack.hh"
+#include "sbnobj/SBND/CRT/CRTSpacePoint.hh"
 #include "sbnobj/Common/Reco/CRUMBSResult.h"
 
 #include "TTree.h"
@@ -70,31 +72,31 @@ namespace sbn {
     void GetMaps(art::Event const& e, std::map<int, int> &trackIDToGenMap, std::map<int, std::string> &genTypeMap,
 		 std::map<int, int> &genCCNCMap, std::map<int, int> &genNuTypeMap);
 
-    art::Ptr<recob::PFParticle> GetSlicePrimary(art::Event const& e, 
-                                                const art::Ptr<recob::Slice> &slice, 
+    art::Ptr<recob::PFParticle> GetSlicePrimary(art::Event const& e,
+                                                const art::Ptr<recob::Slice> &slice,
                                                 const art::ValidHandle<std::vector<recob::Slice> > &handleSlices);
 
-    std::vector<art::Ptr<anab::T0> > GetCRTTrackT0s(art::Event const& e, const art::Ptr<recob::Slice> &slice, 
-                                                    const art::ValidHandle<std::vector<recob::PFParticle> > &handlePFPs,
-                                                    const art::ValidHandle<std::vector<recob::Slice> > &handleSlices);
+    std::vector<anab::T0> GetCRTTrackT0s(art::Event const& e, const art::Ptr<recob::Slice> &slice,
+                                         const art::ValidHandle<std::vector<recob::PFParticle> > &handlePFPs,
+                                         const art::ValidHandle<std::vector<recob::Slice> > &handleSlices);
 
-    std::vector<art::Ptr<anab::T0> > GetCRTHitT0s(art::Event const& e, const art::Ptr<recob::Slice> &slice, 
-                                                  const art::ValidHandle<std::vector<recob::PFParticle> > &handlePFPs,
-                                                  const art::ValidHandle<std::vector<recob::Slice> > &handleSlices);
+    std::vector<anab::T0> GetCRTSPT0s(art::Event const& e, const art::Ptr<recob::Slice> &slice,
+                                      const art::ValidHandle<std::vector<recob::PFParticle> > &handlePFPs,
+                                      const art::ValidHandle<std::vector<recob::Slice> > &handleSlices);
 
-    float GetLongestTrackStoppingChi2Ratio(art::Event const& e, const art::Ptr<recob::Slice> &slice, 
+    float GetLongestTrackStoppingChi2Ratio(art::Event const& e, const art::Ptr<recob::Slice> &slice,
                                            const art::ValidHandle<std::vector<recob::PFParticle> > &handlePFPs,
                                            const art::ValidHandle<std::vector<recob::Slice> > &handleSlices);
 
-    void FillCRTVars(const std::vector<art::Ptr<anab::T0> > &trackT0s, const std::vector<art::Ptr<anab::T0> > &hitT0s);
+    void FillCRTVars(const std::vector<anab::T0> &trackT0s, const std::vector<anab::T0> &hitT0s);
 
     void FillPandoraNuScoreVars(std::map<std::string, float> &propertiesMap);
 
-    std::vector<art::Ptr<recob::Hit> > GetAllSliceHits(art::Event const& e, 
-                                                       const art::Ptr<recob::Slice> &slice, 
+    std::vector<art::Ptr<recob::Hit> > GetAllSliceHits(art::Event const& e,
+                                                       const art::Ptr<recob::Slice> &slice,
                                                        const art::ValidHandle<std::vector<recob::Slice> > &handleSlices);
 
-    void GetTruthMatching(art::Event const& e, const std::vector<art::Ptr<recob::Hit> > &sliceHits, const std::vector<art::Ptr<recob::Hit> > &allHits, 
+    void GetTruthMatching(art::Event const& e, const std::vector<art::Ptr<recob::Hit> > &sliceHits, const std::vector<art::Ptr<recob::Hit> > &allHits,
                           std::map<int, int> &trackIDToGenMap, int &matchedID, double &purity, double &completeness);
 
     int SliceTruthId(std::map<int, float> &purities);
@@ -106,7 +108,7 @@ namespace sbn {
 
     // Module labels
     std::string fMCParticleModuleLabel, fGeneratorModuleLabel, fCosmicModuleLabel, fPFParticleModuleLabel, fHitModuleLabel, fTrackModuleLabel, fSliceModuleLabel, 
-      fFlashMatchModuleLabel, fCRTTrackMatchModuleLabel, fCRTHitMatchModuleLabel, fCalorimetryModuleLabel;
+      fFlashMatchModuleLabel, fCRTTrackMatchModuleLabel, fCRTSPMatchModuleLabel, fCalorimetryModuleLabel;
 
     // MVA location and type for loading
     std::string fMVAName, fMVAFileName, fCCNuMuMVAName, fCCNuMuMVAFileName, fCCNuEMVAName, fCCNuEMVAFileName, fNCMVAName, fNCMVAFileName;
@@ -159,9 +161,9 @@ namespace sbn {
 
     // CRT Track and Hit Matching Variables
     float crt_TrackScore;                // a combination of the DCA and angle between the best matched TPC & CRT tracks
-    float crt_HitScore;                  // the best distance from an extrapolated TPC track to a CRT hit [cm]
+    float crt_SPScore;                   // the best distance from an extrapolated TPC track to a CRT SP [cm]
     float crt_TrackTime;                 // the time associated with the matched CRT track [us]
-    float crt_HitTime;                   // the time associated with the matched CRT hit [us]
+    float crt_SPTime;                    // the time associated with the matched CRT SP [us]
   };
 
 
@@ -180,7 +182,7 @@ namespace sbn {
     fSliceModuleLabel             (p.get<std::string>("SliceModuleLabel")),
     fFlashMatchModuleLabel        (p.get<std::string>("FlashMatchModuleLabel")),
     fCRTTrackMatchModuleLabel     (p.get<std::string>("CRTTrackMatchModuleLabel")),
-    fCRTHitMatchModuleLabel       (p.get<std::string>("CRTHitMatchModuleLabel")),
+    fCRTSPMatchModuleLabel        (p.get<std::string>("CRTSPMatchModuleLabel")),
     fCalorimetryModuleLabel       (p.get<std::string>("CalorimetryModuleLabel")),
     fMVAName                      (p.get<std::string>("MVAName")),
     fMVAFileName                  (p.get<std::string>("MVAFileName")),
@@ -227,9 +229,9 @@ namespace sbn {
           fSliceTree->Branch("pds_FMTime",&pds_FMTime);
 
           fSliceTree->Branch("crt_TrackScore",&crt_TrackScore);
-          fSliceTree->Branch("crt_HitScore",&crt_HitScore);
+          fSliceTree->Branch("crt_SPScore",&crt_SPScore);
           fSliceTree->Branch("crt_TrackTime",&crt_TrackTime);
-          fSliceTree->Branch("crt_HitTime",&crt_HitTime);
+          fSliceTree->Branch("crt_SPTime",&crt_SPTime);
 
           fSliceTree->Branch("eventID",&eventID);
           fSliceTree->Branch("subRunID",&subRunID);
@@ -262,9 +264,9 @@ namespace sbn {
     mvaReader.AddVariable("pds_FMTime",&pds_FMTime);
 
     mvaReader.AddVariable("crt_TrackScore",&crt_TrackScore);
-    mvaReader.AddVariable("crt_HitScore",&crt_HitScore);
+    mvaReader.AddVariable("crt_HitScore",&crt_SPScore);
     mvaReader.AddVariable("crt_TrackTime",&crt_TrackTime);
-    mvaReader.AddVariable("crt_HitTime",&crt_HitTime);
+    mvaReader.AddVariable("crt_HitTime",&crt_SPTime);
 
     cet::search_path searchPath("FW_SEARCH_PATH");
     std::string weightFileFullPath;
@@ -282,7 +284,7 @@ namespace sbn {
 
     pds_FMTotalScore = -999999.; pds_FMPE = -999999.; pds_FMTime = -500.;
 
-    crt_TrackScore = -4.; crt_HitScore = -4.; crt_TrackTime = -3000; crt_HitTime = -3000;
+    crt_TrackScore = -4.; crt_SPScore = -4.; crt_TrackTime = -3000; crt_SPTime = -3000;
 
     slicePDG = 999999;
     matchedType = "";
@@ -396,14 +398,14 @@ namespace sbn {
 
         const std::vector<art::Ptr<larpandoraobj::PFParticleMetadata> > pfpMetaVec = pfpMetadataAssoc.at(primary.key());
         const std::vector<art::Ptr<sbn::SimpleFlashMatch> > pfpFMVec = pfpFMAssoc.at(primary.key());
-        const std::vector<art::Ptr<anab::T0> > sliceCRTTrackT0s = this->GetCRTTrackT0s(e, slice, handlePFPs, handleSlices);
-        const std::vector<art::Ptr<anab::T0> > sliceCRTHitT0s = this->GetCRTHitT0s(e, slice, handlePFPs, handleSlices);
+        const std::vector<anab::T0> sliceCRTTrackT0s = this->GetCRTTrackT0s(e, slice, handlePFPs, handleSlices);
+        const std::vector<anab::T0> sliceCRTSPT0s = this->GetCRTSPT0s(e, slice, handlePFPs, handleSlices);
 
-        this->FillCRTVars(sliceCRTTrackT0s, sliceCRTHitT0s);
+        this->FillCRTVars(sliceCRTTrackT0s, sliceCRTSPT0s);
 
         const art::Ptr<larpandoraobj::PFParticleMetadata> pfpMeta = pfpMetaVec.front();
         std::map<std::string, float> propertiesMap = pfpMeta->GetPropertiesMap();
-      
+
         this->FillPandoraNuScoreVars(propertiesMap);
 
         tpc_StoppingChi2CosmicRatio = this->GetLongestTrackStoppingChi2Ratio(e, slice, handlePFPs, handleSlices);
@@ -426,8 +428,8 @@ namespace sbn {
 	    resultsVec->emplace_back(score, ccnumuscore, ccnuescore, ncscore, bestscore, bestid, tpc_CRFracHitsInLongestTrack, tpc_CRLongestTrackDeflection, 
 				     tpc_CRLongestTrackDirY, std::round(tpc_CRNHitsMax), tpc_NuEigenRatioInSphere, std::round(tpc_NuNFinalStatePfos), 
 				     std::round(tpc_NuNHitsTotal), std::round(tpc_NuNSpacePointsInSphere), tpc_NuVertexY, tpc_NuWeightedDirZ, 
-				     tpc_StoppingChi2CosmicRatio, pds_FMTotalScore, pds_FMPE, pds_FMTime, crt_TrackScore, crt_HitScore, 
-				     crt_TrackTime, crt_HitTime);
+				     tpc_StoppingChi2CosmicRatio, pds_FMTotalScore, pds_FMPE, pds_FMTime, crt_TrackScore, crt_SPScore, 
+				     crt_TrackTime, crt_SPTime);
 	    
 	    util::CreateAssn(*this, e, *resultsVec, slice, *sliceAssns);
 	  }
@@ -456,28 +458,28 @@ namespace sbn {
       }
   }
 
-  void CRUMBS::FillCRTVars(const std::vector<art::Ptr<anab::T0> > &trackT0s, const std::vector<art::Ptr<anab::T0> > &hitT0s)
+  void CRUMBS::FillCRTVars(const std::vector<anab::T0> &trackT0s, const std::vector<anab::T0> &spT0s)
   {
     if (!trackT0s.empty()){
       crt_TrackScore = std::numeric_limits<float>::max();
       for(auto const crttrackmatcht0 : trackT0s)
         {
-          if(crttrackmatcht0->TriggerConfidence() < crt_TrackScore)
+          if(crttrackmatcht0.TriggerConfidence() < crt_TrackScore)
             {
-              crt_TrackScore = crttrackmatcht0->TriggerConfidence();
-              crt_TrackTime = crttrackmatcht0->Time() * 1e-3;
+              crt_TrackScore = crttrackmatcht0.TriggerConfidence();
+              crt_TrackTime  = crttrackmatcht0.Time() * 1e-3;
             }
         }
     }
   
-    if (!hitT0s.empty()){
-      crt_HitScore = std::numeric_limits<float>::max();
-      for(auto const crthitmatcht0 : hitT0s)
+    if (!spT0s.empty()){
+      crt_SPScore = std::numeric_limits<float>::max();
+      for(auto const crtspmatcht0 : spT0s)
         {
-          if(crthitmatcht0->TriggerConfidence() < crt_HitScore)
+          if(crtspmatcht0.TriggerConfidence() < crt_SPScore)
             {
-              crt_HitScore = crthitmatcht0->TriggerConfidence();
-              crt_HitTime = crthitmatcht0->Time() * 1e-3;
+              crt_SPScore = crtspmatcht0.TriggerConfidence();
+              crt_SPTime  = crtspmatcht0.Time() * 1e-3;
             }
         }
     }
@@ -625,17 +627,17 @@ namespace sbn {
       completeness = sliceHitMap[matchedID] / (float) totalTrueHits;
   }
 
-  std::vector<art::Ptr<anab::T0> > CRUMBS::GetCRTTrackT0s(art::Event const& e, const art::Ptr<recob::Slice> &slice, const art::ValidHandle<std::vector<recob::PFParticle> > &handlePFPs,
-                                                          const art::ValidHandle<std::vector<recob::Slice> > &handleSlices)
+  std::vector<anab::T0> CRUMBS::GetCRTTrackT0s(art::Event const& e, const art::Ptr<recob::Slice> &slice, const art::ValidHandle<std::vector<recob::PFParticle> > &handlePFPs,
+                                               const art::ValidHandle<std::vector<recob::Slice> > &handleSlices)
   {
-    std::vector<art::Ptr<anab::T0> > t0Vec;
+    std::vector<anab::T0> t0Vec;
 
     art::Handle<std::vector<recob::Track> > handleTracks;
     e.getByLabel(fTrackModuleLabel, handleTracks);
 
     art::FindManyP<recob::PFParticle> slicePFPAssn(handleSlices,e,fSliceModuleLabel);
     art::FindManyP<recob::Track> pfpTrackAssn(handlePFPs,e,fTrackModuleLabel);
-    art::FindManyP<anab::T0> trackT0Assn(handleTracks,e,fCRTTrackMatchModuleLabel);
+    art::FindOneP<sbnd::crt::CRTTrack, anab::T0> trackT0Assn(handleTracks,e,fCRTTrackMatchModuleLabel);
 
     const std::vector<art::Ptr<recob::PFParticle> > pfps = slicePFPAssn.at(slice.key());
   
@@ -651,24 +653,29 @@ namespace sbn {
 
         const art::Ptr<recob::Track> track = tracks.front();
 
-        const std::vector<art::Ptr<anab::T0> > t0s = trackT0Assn.at(track.key());
-        t0Vec.insert(t0Vec.end(), t0s.begin(), t0s.end());
+	const art::Ptr<sbnd::crt::CRTTrack> crttrack = trackT0Assn.at(track.key());
+
+	if(crttrack.isNonnull())
+	  {
+	    const anab::T0 t0 = trackT0Assn.data(track.key()).ref();
+	    t0Vec.push_back(t0);
+	  }
       }
   
     return t0Vec;
   }
 
-  std::vector<art::Ptr<anab::T0> > CRUMBS::GetCRTHitT0s(art::Event const& e, const art::Ptr<recob::Slice> &slice, const art::ValidHandle<std::vector<recob::PFParticle> > &handlePFPs,
-                                                        const art::ValidHandle<std::vector<recob::Slice> > &handleSlices)
+  std::vector<anab::T0> CRUMBS::GetCRTSPT0s(art::Event const& e, const art::Ptr<recob::Slice> &slice, const art::ValidHandle<std::vector<recob::PFParticle> > &handlePFPs,
+                                            const art::ValidHandle<std::vector<recob::Slice> > &handleSlices)
   {
-    std::vector<art::Ptr<anab::T0> > t0Vec;
+    std::vector<anab::T0> t0Vec;
 
     art::Handle<std::vector<recob::Track> > handleTracks;
     e.getByLabel(fTrackModuleLabel, handleTracks);
 
     art::FindManyP<recob::PFParticle> slicePFPAssn(handleSlices,e,fSliceModuleLabel);
     art::FindManyP<recob::Track> pfpTrackAssn(handlePFPs,e,fTrackModuleLabel);
-    art::FindManyP<anab::T0> trackT0Assn(handleTracks,e,fCRTHitMatchModuleLabel);
+    art::FindOneP<sbnd::crt::CRTSpacePoint, anab::T0> trackT0Assn(handleTracks,e,fCRTSPMatchModuleLabel);
 
     const std::vector<art::Ptr<recob::PFParticle> > pfps = slicePFPAssn.at(slice.key());
   
@@ -684,8 +691,13 @@ namespace sbn {
 
         const art::Ptr<recob::Track> track = tracks.front();
 
-        const std::vector<art::Ptr<anab::T0> > t0s = trackT0Assn.at(track.key());
-        t0Vec.insert(t0Vec.end(), t0s.begin(), t0s.end());
+	const art::Ptr<sbnd::crt::CRTSpacePoint> crtsp = trackT0Assn.at(track.key());
+
+	if(crtsp.isNonnull())
+	  {
+	    const anab::T0 t0 = trackT0Assn.data(track.key()).ref();
+	    t0Vec.push_back(t0);
+	  }
       }
   
     return t0Vec;
diff --git a/sbncode/CosmicID/sbn_crumbs_producer.fcl b/sbncode/CosmicID/sbn_crumbs_producer.fcl
index a71af130d..08600df57 100644
--- a/sbncode/CosmicID/sbn_crumbs_producer.fcl
+++ b/sbncode/CosmicID/sbn_crumbs_producer.fcl
@@ -13,8 +13,8 @@ crumbs_sbnd:
   TrackModuleLabel:			"pandoraTrack"
   SliceModuleLabel:			"pandora"
   FlashMatchModuleLabel:		"fmatch"
-  CRTTrackMatchModuleLabel:		"crttrackt0"
-  CRTHitMatchModuleLabel:		"crthitt0"
+  CRTTrackMatchModuleLabel:		"crttrackmatching"
+  CRTSPMatchModuleLabel:		"crtspacepointmatching"
   CalorimetryModuleLabel:		"pandoraCalo"
   
   MVAName:                              "BDT"
diff --git a/sbncode/CosmicID/training_scripts/CRUMBSTMVADriver.C b/sbncode/CosmicID/training_scripts/CRUMBSTMVADriver.C
index 6cc90416f..bc757d4a4 100644
--- a/sbncode/CosmicID/training_scripts/CRUMBSTMVADriver.C
+++ b/sbncode/CosmicID/training_scripts/CRUMBSTMVADriver.C
@@ -91,9 +91,9 @@ void TrainCRUMBSInstance(const TString outDirName, TTree *inputTree,
   dataloader->AddVariable("pds_FMTime","FM Time","#mu s",'F');
 
   dataloader->AddVariable("crt_TrackScore","CRT Track Match Score","",'F');
-  dataloader->AddVariable("crt_HitScore","CRT Hit Match Score","",'F');
+  dataloader->AddVariable("crt_SPScore","CRT SpacePoint Match Score","",'F');
   dataloader->AddVariable("crt_TrackTime","CRT Track Match Time","#mu s",'F');
-  dataloader->AddVariable("crt_HitTime","CRT Hit Match Time","#mu s",'F');
+  dataloader->AddVariable("crt_SPTime","CRT SpacePoint Match Time","#mu s",'F');
 
   dataloader->AddSignalTree    (inputTree, signalWeight);
   dataloader->AddBackgroundTree(inputTree, backgroundWeight);
diff --git a/sbncode/TPCReco/CMakeLists.txt b/sbncode/TPCReco/CMakeLists.txt
index 0d2338bd6..9ff783574 100644
--- a/sbncode/TPCReco/CMakeLists.txt
+++ b/sbncode/TPCReco/CMakeLists.txt
@@ -2,6 +2,9 @@ add_subdirectory(TrackHit)
 add_subdirectory(TrackSplit)
 add_subdirectory(VertexStub)
 add_subdirectory(CalorimetryAnalysis)
+if (TensorFlow_FOUND)
+  add_subdirectory(CNNHitClassification)
+endif()
 
 
 cet_build_plugin(NuVertexChargeTree art::module
diff --git a/sbncode/TPCReco/CNNHitClassification/CMakeLists.txt b/sbncode/TPCReco/CNNHitClassification/CMakeLists.txt
new file mode 100644
index 000000000..709fdc440
--- /dev/null
+++ b/sbncode/TPCReco/CNNHitClassification/CMakeLists.txt
@@ -0,0 +1,45 @@
+cet_build_plugin(CNNID art::EDProducer
+  LIBRARIES PRIVATE
+    larrecodnn::ImagePatternAlgs_Tensorflow_PointIdAlg
+    lardata::ArtDataHelper
+    lardata::DetectorClocksService
+    lardata::DetectorPropertiesService
+    lardata::AssociationUtil
+    lardataobj::RecoBase
+    sbnobj::Common_Reco
+    larcoreobj::SimpleTypesAndConstants
+    art::Framework_Services_System_TriggerNamesService_service
+    art::Framework_Services_Registry
+    canvas::canvas
+    messagefacility::MF_MessageLogger
+    fhiclcpp::types
+    fhiclcpp::fhiclcpp
+    cetlib_except::cetlib_except
+)
+
+# cet_build_plugin(CheckCNNScoreAllHits art::EDProducer
+#   LIBRARIES PRIVATE
+#     larrecodnn::ImagePatternAlgs_Tensorflow_PointIdAlg
+#     lardata::ArtDataHelper
+#     lardata::DetectorClocksService
+#     lardata::DetectorPropertiesService
+#     lardata::AssociationUtil
+#     lardataobj::RecoBase
+#     sbnobj::Common_Reco
+#     larcoreobj::SimpleTypesAndConstants
+#     art::Framework_Services_System_TriggerNamesService_service
+#     art::Framework_Services_Registry
+#     canvas::canvas
+#     messagefacility::MF_MessageLogger
+#     fhiclcpp::types
+#     fhiclcpp::fhiclcpp
+#     cetlib_except::cetlib_except
+#     art_root_io::TFileService_service
+#     art_root_io::tfile_support
+#     art_root_io::art_root_io
+#     art_root_io::dict
+# )
+
+install_headers()
+install_fhicl()
+install_source()
\ No newline at end of file
diff --git a/sbncode/TPCReco/CNNHitClassification/CNNID_module.cc b/sbncode/TPCReco/CNNHitClassification/CNNID_module.cc
new file mode 100644
index 000000000..d85281f68
--- /dev/null
+++ b/sbncode/TPCReco/CNNHitClassification/CNNID_module.cc
@@ -0,0 +1,406 @@
+/////////////////////////////////////////////////////////////////////////////////////
+// Class:       CNNID
+// Module Type: producer
+// File:        CNNID_module.cc
+//
+// apply CNN to hits associated with cluster to get track/shower/noise/endMichel score
+// score of a PFP is the average score of all associated hits
+// skip clear cosmic PFPs
+//
+// M.Jung munjung@uchicago.edu
+/////////////////////////////////////////////////////////////////////////////////////
+
+#include "larrecodnn/ImagePatternAlgs/Tensorflow/PointIdAlg/PointIdAlg.h"
+#include "lardata/ArtDataHelper/MVAWriter.h"
+#include "lardata/DetectorInfoServices/DetectorClocksService.h"
+#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
+#include "lardata/Utilities/AssociationUtil.h"
+#include "lardataobj/RecoBase/Cluster.h"
+#include "lardataobj/RecoBase/Hit.h"
+#include "lardataobj/RecoBase/Track.h"
+#include "lardataobj/RecoBase/Wire.h"
+#include "lardataobj/RecoBase/PFParticle.h"
+#include "lardataobj/RecoBase/PFParticleMetadata.h"
+#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
+
+#include "art/Framework/Core/EDProducer.h"
+#include "art/Framework/Core/ModuleMacros.h"
+#include "art/Framework/Principal/Event.h"
+#include "art/Framework/Principal/Handle.h"
+#include "art/Framework/Principal/Run.h"
+#include "art/Framework/Services/Registry/ServiceHandle.h"
+#include "art/Framework/Services/System/TriggerNamesService.h"
+#include "canvas/Persistency/Common/Assns.h"
+#include "canvas/Persistency/Common/Ptr.h"
+#include "canvas/Persistency/Common/FindManyP.h"
+#include "canvas/Utilities/InputTag.h"
+#include "messagefacility/MessageLogger/MessageLogger.h"
+#include "fhiclcpp/types/Atom.h"
+#include "fhiclcpp/types/Comment.h"
+#include "fhiclcpp/types/Name.h"
+#include "fhiclcpp/types/Sequence.h"
+#include "fhiclcpp/types/Table.h"
+#include "fhiclcpp/ParameterSet.h"
+#include "cetlib_except/exception.h"
+#include "sbnobj/Common/Reco/CNNScore.h"
+
+#include <iostream>
+#include <memory>
+#include <string>
+#include <unordered_map>
+#include <utility>
+#include <vector>
+
+namespace sbn {
+
+  class CNNID : public art::EDProducer {
+  public:
+
+    struct Config {
+      using Name = fhicl::Name;
+      using Comment = fhicl::Comment;
+
+      fhicl::Table<nnet::PointIdAlg::Config> PointIdAlg{Name("PointIdAlg")};
+      fhicl::Atom<size_t> BatchSize{Name("BatchSize"),
+                                    Comment("number of samples processed in one batch")};
+
+      fhicl::Atom<art::InputTag> WireLabel{
+        Name("WireLabel"),
+        Comment("tag of deconvoluted ADC on wires (recob::Wire)")};
+
+      fhicl::Atom<art::InputTag> HitModuleLabel{
+        Name("HitModuleLabel"),
+        Comment("tag of hits to be EM/track / Michel tagged")};
+
+      fhicl::Atom<art::InputTag> ClusterModuleLabel{
+        Name("ClusterModuleLabel"),
+        Comment("tag of clusters")};
+
+      fhicl::Atom<art::InputTag> PFParticleModuleLabel{
+        Name("PFParticleModuleLabel"),
+        Comment("tag of PFParticles")};
+
+      fhicl::Atom<bool> SkipClearCosmics{
+        Name("SkipClearCosmics"),
+        Comment("skip clear cosmic PFPs")};
+
+      fhicl::Atom<bool> DoMichel{
+        Name("DoMichel"),
+        Comment("get Michel Scores")};
+    
+      fhicl::Sequence<int> MichelRegionSize{
+        Name("MichelRegionSize"),
+        Comment("size of region around cluster end to get Michel Scores")};
+
+      fhicl::Atom<bool> DoPFP{
+        Name("DoPFP"),
+        Comment("get PFP Scores")};
+
+      fhicl::Atom<int> SparseLengthCut{
+        Name("SparseLengthCut"),
+        Comment("sparsely apply CNN if the number of hits is greater than this number")};
+
+      fhicl::Atom<int> SparseRate{
+        Name("SparseRate"),
+        Comment("rate to sparsely apply CNN")};
+
+      fhicl::Sequence<int> Views{
+        Name("Views"),
+        Comment("tag clusters in selected views only, or in all views if empty list")};
+    };
+    using Parameters = art::EDProducer::Table<Config>;
+    explicit CNNID(Parameters const& p);
+
+    CNNID(CNNID const&) = delete;
+    CNNID(CNNID&&) = delete;
+    CNNID& operator=(CNNID const&) = delete;
+    CNNID& operator=(CNNID&&) = delete;
+
+  private:
+    void produce(art::Event& e) override;
+    float getAvgScore(std::vector<float> scores);
+    size_t fBatchSize;
+    nnet::PointIdAlg fPointIdAlg;
+    anab::MVAWriter<4> fMVAWriter; 
+    art::InputTag fWireProducerLabel;
+    art::InputTag fHitModuleLabel;
+    art::InputTag fClusterModuleLabel;
+    art::InputTag fPFParticleModuleLabel;
+    bool fSkipClearCosmics;
+    bool fDoMichel;
+    std::vector<int> fMichelRegionSize;
+    bool fDoPFP;
+    int fSparseLengthCut;
+    int fSparseRate;
+    std::vector<int> fViews;
+
+  };
+  // ------------------------------------------------------
+
+  CNNID::CNNID(CNNID::Parameters const& config)
+    : EDProducer{config}
+    , fBatchSize(config().BatchSize())
+    , fPointIdAlg(config().PointIdAlg())
+    , fMVAWriter(producesCollector(), "cnnid")
+    , fWireProducerLabel(config().WireLabel())
+    , fHitModuleLabel(config().HitModuleLabel())
+    , fClusterModuleLabel(config().ClusterModuleLabel())
+    , fPFParticleModuleLabel(config().PFParticleModuleLabel())
+    , fSkipClearCosmics(config().SkipClearCosmics())
+    , fDoMichel(config().DoMichel())
+    , fMichelRegionSize(config().MichelRegionSize())
+    , fDoPFP(config().DoPFP())
+    , fSparseLengthCut(config().SparseLengthCut())
+    , fSparseRate(config().SparseRate())
+    , fViews(config().Views())
+  {
+    produces<std::vector<PFPCNNScore>>();
+    produces<art::Assns<recob::PFParticle, PFPCNNScore>>();
+    fMVAWriter.produces_using<recob::Hit>();
+  }
+  // ------------------------------------------------------
+
+  void
+  CNNID::produce(art::Event& evt)
+  {
+    mf::LogVerbatim("CNNID") << "next event: " << evt.run() << " / " << evt.id().event();
+
+    auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
+    auto const detProp =
+        art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);
+
+    auto wireHandle = evt.getValidHandle<std::vector<recob::Wire>>(fWireProducerLabel);
+
+    auto hitListHandle = evt.getValidHandle<std::vector<recob::Hit>>(fHitModuleLabel);
+    std::vector<art::Ptr<recob::Hit>> hitPtrList;
+    art::fill_ptr_vector(hitPtrList, hitListHandle);
+
+    auto pfpListHandle = evt.getValidHandle<std::vector<recob::PFParticle>>(fPFParticleModuleLabel);
+		std::vector<art::Ptr<recob::PFParticle>> pfpList;
+    art::fill_ptr_vector(pfpList, pfpListHandle);
+
+    auto metaHandle = evt.getValidHandle<std::vector<larpandoraobj::PFParticleMetadata>>(fPFParticleModuleLabel);
+    std::vector<art::Ptr<larpandoraobj::PFParticleMetadata>> metaPtrList;
+    art::fill_ptr_vector(metaPtrList, metaHandle);
+
+    auto cluListHandle = evt.getValidHandle<std::vector<recob::Cluster>>(fClusterModuleLabel);
+    std::vector<art::Ptr<recob::Cluster>> cluPtrList;
+    art::fill_ptr_vector(cluPtrList, cluListHandle);
+
+    art::FindManyP<larpandoraobj::PFParticleMetadata> metaFMpfp(pfpListHandle, evt, fPFParticleModuleLabel);
+    art::FindManyP<recob::Cluster> cluFMpfp(pfpListHandle, evt, fPFParticleModuleLabel);
+    art::FindManyP<recob::PFParticle> pfpFMclu(cluListHandle, evt, fClusterModuleLabel);
+    art::FindManyP<recob::Hit> hitFMclu(cluListHandle, evt, fClusterModuleLabel);
+    art::FindManyP<recob::Hit> hitFMpfp(pfpListHandle, evt, fPFParticleModuleLabel);
+    art::FindManyP<recob::Cluster> cluFMhit(hitListHandle, evt, fClusterModuleLabel);
+
+    auto hitID = fMVAWriter.initOutputs<recob::Hit>(
+      fHitModuleLabel, hitPtrList.size(), fPointIdAlg.outputLabels());
+
+    auto pfpScores = std::make_unique<std::vector<PFPCNNScore>>();
+    auto pfpAssns = std::make_unique<art::Assns<recob::PFParticle, PFPCNNScore>>();
+
+    // loop over PFPs
+    for (const art::Ptr<recob::PFParticle> &pfp : pfpList){
+      std::vector<float> pfpTrackScore; 
+      std::vector<float> pfpShowerScore; 
+      std::vector<float> pfpNoiseScore;
+      std::vector<float> pfpMichelScore;
+      std::vector<float> pfpEndMichelScore;
+      float pfpAvgTrackScore;
+      float pfpAvgShowerScore;
+      float pfpAvgNoiseScore;
+      float pfpAvgMichelScore;
+      float pfpAvgEndMichelScore;
+      int nClusters = 0; 
+
+      if (fSkipClearCosmics) { // skip clear cosmic PFPs
+        std::vector<art::Ptr<larpandoraobj::PFParticleMetadata>> metas = metaFMpfp.at(pfp.key());
+        auto const &properties = metas[0]->GetPropertiesMap();
+        if (properties.count("IsClearCosmic")) {
+        //    continue;
+        // TODO: assign nan values instead of skipping?
+            pfpAvgTrackScore = std::numeric_limits<float>::signaling_NaN();
+            pfpAvgShowerScore = std::numeric_limits<float>::signaling_NaN();
+            pfpAvgNoiseScore = std::numeric_limits<float>::signaling_NaN();
+            pfpAvgMichelScore = std::numeric_limits<float>::signaling_NaN();
+            pfpAvgEndMichelScore = std::numeric_limits<float>::signaling_NaN();
+            pfpScores->emplace_back(pfpAvgTrackScore, pfpAvgShowerScore, pfpAvgNoiseScore, pfpAvgMichelScore, pfpAvgEndMichelScore, nClusters); //
+            util::CreateAssn(*this, evt, *pfpScores, pfp, *pfpAssns);
+            continue;
+        }
+      }
+
+      // loop over clusters
+      auto const &clusters = cluFMpfp.at(pfp->Self());
+      for (auto const &cluster : clusters){
+        nClusters++;
+        unsigned int cluView = cluster->Plane().Plane;
+        unsigned int cluCryo = cluster->Plane().Cryostat;
+        unsigned int cluTPC = cluster->Plane().TPC;
+        fPointIdAlg.setWireDriftData(clockData, detProp, *wireHandle, cluView, cluTPC, cluCryo);
+
+        float cluTrackScore = 0.;
+        float cluShowerScore = 0.;
+        float cluNoiseScore = 0.;
+        float cluMichelScore = 0.;
+        float cluEndMichelScore = 0.;
+
+        if (fDoMichel) { // Michel scores for hits around end of cluster
+          int clueEndWire = cluster->EndWire();
+          auto clueEndTick = cluster->EndTick();
+          std::vector<art::Ptr<recob::Hit>> michelPtrList;
+          for (auto const& hit : hitPtrList) {
+            unsigned int hitView = hit->WireID().Plane;
+            int hitWireID = hit->WireID().Wire;
+            auto hitPeakTime = hit->PeakTime();
+
+            if (hitView != cluView) continue; // same plane
+            if ((std::abs(hitWireID-clueEndWire) > fMichelRegionSize[0]) || (std::abs(hitPeakTime-clueEndTick) > fMichelRegionSize[1])) continue;
+            bool isValidMichelHit = true;
+            auto const &michelClusters = cluFMhit.at(hit.key());
+            for (auto const &michelCluster : michelClusters){
+              // exclude if hit is on the parent cluster 
+              if (michelCluster.key() == cluster.key()) {
+                isValidMichelHit = false;
+                break;
+              }
+              // exclude if hit is on a cosmic cluster
+              auto const &michelPFPs = pfpFMclu.at(michelCluster.key());
+              for (auto const &michelPFP : michelPFPs){
+                std::vector<art::Ptr<larpandoraobj::PFParticleMetadata>> michelMetas = metaFMpfp.at(michelPFP.key());
+                auto const &michelProperties = michelMetas[0]->GetPropertiesMap();
+                if (michelProperties.count("IsClearCosmic")) {
+                  isValidMichelHit = false;
+                  break;
+                }
+              }
+            }
+            if (isValidMichelHit) michelPtrList.push_back(hit);
+          }
+          std::vector<std::pair<unsigned int, float>> michelHits;
+          std::vector<size_t> keys;
+          for (auto const& hit : michelPtrList){
+            michelHits.emplace_back(hit->WireID().Wire, hit->PeakTime());
+            keys.push_back(hit.key());
+          }
+          auto batchOut = fPointIdAlg.predictIdVectors(michelHits);
+          if (michelHits.size() != batchOut.size()) {
+              throw cet::exception("CNNID") << "Size mismatch between input and output vectors";
+          }
+          size_t nMichelHits = michelHits.size();
+          for (size_t h = 0; h < nMichelHits; h++) {
+            fMVAWriter.setOutput(hitID, keys[h], batchOut[h]);
+            cluEndMichelScore += batchOut[h][3];
+          }
+          if (nMichelHits != 0) cluEndMichelScore = cluEndMichelScore/nMichelHits;
+        }
+        else cluEndMichelScore = std::numeric_limits<float>::signaling_NaN();
+        // Michel score end
+
+        // non clear cosmic pfp track/shower/noise score
+        if (fDoPFP) {
+          auto &cluHitPtrList = hitFMclu.at(cluster.key());
+          std::vector<art::Ptr<recob::Hit>> cluHitPtrListApply; 
+          int nCluHits = cluHitPtrList.size();
+          if (nCluHits > fSparseLengthCut) {
+            for (int i = 0; i < nCluHits; i++){
+              if (i % fSparseRate != 0) continue;
+              cluHitPtrListApply.push_back(cluHitPtrList[i]);
+            }
+          }
+          else cluHitPtrListApply = cluHitPtrList;
+          std::vector<std::pair<unsigned int, float>> cluApplyHits;
+          std::vector<size_t> keys;
+          for (auto const& hit : cluHitPtrListApply){
+            cluApplyHits.emplace_back(hit->WireID().Wire, hit->PeakTime());
+            keys.push_back(hit.key());
+          }
+          size_t nCluApplyHits = cluApplyHits.size();
+          auto batchOut = fPointIdAlg.predictIdVectors(cluApplyHits);
+          if (nCluApplyHits != batchOut.size()) {
+            throw cet::exception("CNNID") << "Size mismatch between input and output vectors";
+          }
+          for (size_t h = 0; h < nCluApplyHits; h++) {
+            fMVAWriter.setOutput(hitID, keys[h], batchOut[h]);
+            cluTrackScore += batchOut[h][0];
+            cluShowerScore += batchOut[h][1];
+            cluNoiseScore += batchOut[h][2];
+            cluMichelScore += batchOut[h][3];
+          }
+          if (nCluApplyHits !=0) {
+            cluTrackScore = cluTrackScore/nCluApplyHits;
+            cluShowerScore = cluShowerScore/nCluApplyHits;
+            cluNoiseScore = cluNoiseScore/nCluApplyHits;
+            cluMichelScore = cluMichelScore/nCluApplyHits;
+          }
+        }
+        else {
+          cluTrackScore = std::numeric_limits<float>::signaling_NaN();
+          cluShowerScore = std::numeric_limits<float>::signaling_NaN();
+          cluNoiseScore = std::numeric_limits<float>::signaling_NaN();
+          cluMichelScore = std::numeric_limits<float>::signaling_NaN();
+        }
+        // cluster scores end
+
+        pfpTrackScore.push_back(cluTrackScore);
+        pfpShowerScore.push_back(cluShowerScore);
+        pfpNoiseScore.push_back(cluNoiseScore);
+        pfpMichelScore.push_back(cluMichelScore);
+        pfpEndMichelScore.push_back(cluEndMichelScore);
+      }
+
+      pfpAvgTrackScore = getAvgScore(pfpTrackScore);
+      pfpAvgShowerScore = getAvgScore(pfpShowerScore);
+      pfpAvgNoiseScore = getAvgScore(pfpNoiseScore);
+      pfpAvgMichelScore = getAvgScore(pfpMichelScore);
+      pfpAvgEndMichelScore = getAvgScore(pfpEndMichelScore);
+
+      // std::cout << "PFP ID: " << pfp->Self() << std::endl;
+      // std::cout << " average shower score: " << pfpAvgShowerScore << std::endl;
+      // std::cout << " average noise score: " << pfpAvgNoiseScore << std::endl;
+      // std::cout << " average michel score: " << pfpAvgMichelScore << std::endl;
+      // std::cout << " average end michel score: " << pfpAvgEndMichelScore << std::endl;
+      // std::cout << "------------------------------------------------------" << std::endl;
+      pfpScores->emplace_back(pfpAvgTrackScore, pfpAvgShowerScore, pfpAvgNoiseScore, pfpAvgMichelScore, pfpAvgEndMichelScore, nClusters); //
+      util::CreateAssn(*this, evt, *pfpScores, pfp, *pfpAssns);
+    }
+    
+    fMVAWriter.saveOutputs(evt);
+    evt.put(std::move(pfpScores));
+    evt.put(std::move(pfpAssns));
+  }
+  // ------------------------------------------------------
+  // get average score, drop outliers
+  float CNNID::getAvgScore(std::vector<float> scores){
+    // good score means a positive score (default score is set to 0)
+    // if no good score, return NaN
+    // if only 1 good score, use that score as "average"
+    // if 2 good scores, use the average of the 2 scores
+    // if 3 or more good scores, remove outliers if any and use the average of the remaining scores
+    int nScores = scores.size();
+    if (nScores == 0) return std::numeric_limits<float>::signaling_NaN();
+    int nGoodScore = 0;
+    std::vector<float> goodScores;
+    for (int i = 0; i < nScores; i++) {
+      if (scores[i] == 0.) continue;
+      goodScores.push_back(scores[i]);
+      nGoodScore++;
+    }
+    if (nGoodScore == 0) return std::numeric_limits<float>::signaling_NaN();
+    else if (nGoodScore == 1) return goodScores[0];
+    else if (nGoodScore == 2) return (goodScores[0]+goodScores[1])/2.;
+    else {
+      std::sort(goodScores.begin(), goodScores.end());
+      float diffLow = goodScores[1] - goodScores[0];
+      float diffHigh = goodScores[goodScores.size()-1] - goodScores[goodScores.size()-2];
+      if (diffLow > 0.5) return (goodScores[1]+goodScores[2])/2.; // remove lowest score (outlier)
+      else if (diffHigh > 0.5) return (goodScores[0]+goodScores[1])/2.; // remove highest score (outlier)
+      else return (goodScores[0]+goodScores[1]+goodScores[2])/3.;
+    }
+  }
+
+  // ------------------------------------------------------
+
+  DEFINE_ART_MODULE(CNNID)
+}
diff --git a/sbncode/TPCReco/CNNHitClassification/pointidalg_sbnd.fcl b/sbncode/TPCReco/CNNHitClassification/pointidalg_sbnd.fcl
new file mode 100644
index 000000000..2a6957687
--- /dev/null
+++ b/sbncode/TPCReco/CNNHitClassification/pointidalg_sbnd.fcl
@@ -0,0 +1,21 @@
+#include "dataprovider.fcl"
+
+BEGIN_PROLOG
+
+pointidalg_sbnd: @local::standard_dataprovideralg 
+pointidalg_sbnd.NNetModelFile:     "/cvmfs/sbnd.opensciencegrid.org/products/sbnd/sbnd_data/v01_22_00/CNNHitClassification/CNNID_model.pb"
+pointidalg_sbnd.NNetOutputs:       ["track", "em", "none", "michel"]
+pointidalg_sbnd.CalorimetryAlg:    @local::sbnd_calorimetryalgmc
+pointidalg_sbnd.CalibrateAmpl:     false
+pointidalg_sbnd.CalibrateLifetime: false
+pointidalg_sbnd.PatchSizeW:        60
+pointidalg_sbnd.PatchSizeD:        60
+pointidalg_sbnd.DriftWindow:       4
+pointidalg_sbnd.DownscaleFn:       "mean"
+pointidalg_sbnd.DownscaleFullView: true
+pointidalg_sbnd.AdcMin: 0 
+pointidalg_sbnd.AdcMax: 200
+pointidalg_sbnd.OutMin: 0
+pointidalg_sbnd.OutMax: 1
+
+END_PROLOG
\ No newline at end of file
diff --git a/sbncode/TPCReco/CNNHitClassification/run_cnnid_sbnd.fcl b/sbncode/TPCReco/CNNHitClassification/run_cnnid_sbnd.fcl
new file mode 100644
index 000000000..8788cfb2b
--- /dev/null
+++ b/sbncode/TPCReco/CNNHitClassification/run_cnnid_sbnd.fcl
@@ -0,0 +1,64 @@
+#include "services_sbnd.fcl"
+#include "simulationservices_sbnd.fcl"
+#include "calorimetry_sbnd.fcl"
+#include "pointidalg_sbnd.fcl"
+
+process_name: CNNID
+
+services:
+{
+  TFileService: { fileName: "reco_hist.root" }
+  MemoryTracker:     {}
+  TimeTracker:       {}
+  RandomNumberGenerator: {} 
+  message:              @local::sbnd_message_services_prod_debug
+  FileCatalogMetadata:  @local::art_file_catalog_mc
+             			@table::sbnd_simulation_services
+}
+
+source:
+{
+  module_type: RootInput
+  maxEvents:  -1
+}
+
+
+physics:
+{
+  producers:
+  {
+    cnnid: {
+      module_type: CNNID
+      BatchSize:                    256 
+      WireLabel:                    "caldata"
+      HitModuleLabel:               "gaushit"
+      ClusterModuleLabel:           "pandora"
+      PFParticleModuleLabel:        "pandora"
+      SkipClearCosmics:             true
+      DoMichel:                     true
+      MichelRegionSize:             [30, 120]
+      SparseLengthCut:              150
+      DoPFP:                        true
+      SparseRate:                   5
+      Views:                        []
+      PointIdAlg: @local::pointidalg_sbnd
+    }
+  }
+
+  reco: [ cnnid ]
+  stream1:  [ out1 ]
+  trigger_paths: [reco]
+  end_paths:     [stream1]
+}
+
+outputs:
+{
+ out1:
+ {
+   module_type: RootOutput
+  #  fileName:    "/pnfs/sbnd/scratch/users/munjung/cnnid_job/%ifb_cnn.root"
+   fileName:    "%ifb_cnn.root"
+   dataTier:    "full-reconstructed"
+   fastCloning: true
+ }
+}
diff --git a/sbncode/TPCReco/CNNHitClassification/sbn_cnnid.fcl b/sbncode/TPCReco/CNNHitClassification/sbn_cnnid.fcl
new file mode 100644
index 000000000..784e71b28
--- /dev/null
+++ b/sbncode/TPCReco/CNNHitClassification/sbn_cnnid.fcl
@@ -0,0 +1,23 @@
+#include "pointidalg_sbnd.fcl"
+#include "channelstatus_sbnd.fcl"
+
+BEGIN_PROLOG
+
+cnnid_sbnd: {
+  module_type: CNNID
+  BatchSize:                    256 
+  WireLabel:                    "caldata"
+  HitModuleLabel:               "gaushit"
+  ClusterModuleLabel:           "pandora"
+  PFParticleModuleLabel:        "pandora"
+  SkipClearCosmics:             true
+  DoMichel:                     true
+  MichelRegionSize:             [30, 120]
+  SparseLengthCut:              150
+  DoPFP:                        true
+  SparseRate:                   5
+  Views:                        []
+  PointIdAlg: @local::pointidalg_sbnd
+}
+
+END_PROLOG