diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5bddf043e..f751a7de7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -69,6 +69,7 @@ find_package(lardata REQUIRED )
 find_package(larsim REQUIRED )
 find_package(larevt REQUIRED )
 find_package(larreco REQUIRED )
+find_package(larrecodnn REQUIRED )
 find_package(larana REQUIRED )
 find_package(larpandora REQUIRED )
 find_package(larpandoracontent REQUIRED )
diff --git a/fcl/caf/cafmaker_defs.fcl b/fcl/caf/cafmaker_defs.fcl
index 50de910f3..a71f8cbad 100644
--- a/fcl/caf/cafmaker_defs.fcl
+++ b/fcl/caf/cafmaker_defs.fcl
@@ -319,6 +319,9 @@ cafmaker.SystWeightLabels: ["genieweight", "fluxweight"]
 cafmaker.SaveGENIEEventRecord: true # save GENIE event record by default. Turn this off for data cafmaker fcl
 cafmaker.TPCPMTBarycenterMatchLabel: "tpcpmtbarycentermatch"
 cafmaker.TrackHitFillRREndCut: 30 # include entire PID region
+cafmaker.NuGraphSliceHitLabel: "nuslhits"
+cafmaker.NuGraphFilterLabel: "NuGraph:filter"
+cafmaker.NuGraphSemanticLabel: "NuGraph:semantic"
 
 # Add CAFMaker to the list of producers
 caf_preprocess_producers.cafmaker: @local::cafmaker
diff --git a/fcl/reco/Definitions/stage1_icarus_defs.fcl b/fcl/reco/Definitions/stage1_icarus_defs.fcl
index 22b46744c..2906a4e4e 100644
--- a/fcl/reco/Definitions/stage1_icarus_defs.fcl
+++ b/fcl/reco/Definitions/stage1_icarus_defs.fcl
@@ -22,6 +22,7 @@
 #include "supera_modules.fcl"
 #include "crtpmtmatching_parameters.fcl"
 #include "tpcpmtbarycentermatch_config.fcl"
+#include "nugraph_icarus.fcl"
 
 BEGIN_PROLOG
 
@@ -80,6 +81,13 @@ icarus_stage1_producers:
 
   tpcpmtbarycentermatchCryoE: @local::data_tpcpmtbarycentermatchproducer_east
   tpcpmtbarycentermatchCryoW: @local::data_tpcpmtbarycentermatchproducer_west
+
+  ## NuGraph
+  nuslhitsCryoE: @local::nuslhitsCryoE
+  nuslhitsCryoW: @local::nuslhitsCryoW
+  NuGraphCryoE: @local::NuGraphCryoE
+  NuGraphCryoW: @local::NuGraphCryoW
+
 }
 
 icarus_stage1_filters:
@@ -224,13 +232,19 @@ icarus_reco_cluster3DCryoW:        [ cluster3DCryoW ]
 
 icarus_reco_cluster3DCryoE:        [ cluster3DCryoE ]
 
-icarus_reco_pandoraGausCryoW:      [ pandoraGausCryoW,
+icarus_reco_pandoraGausCryoW:      [ pandoraGausCryoW ]
+
+icarus_reco_pandoraGausCryoE:      [ pandoraGausCryoE ]
+
+icarus_reco_postPandoraGausCryoW:  [ nuslhitsCryoW,
+                                     NuGraphCryoW,
                                      pandoraTrackGausCryoW,
                                      pandoraKalmanTrackGausCryoW,
                                      SBNShowerGausCryoW
                                    ]
 
-icarus_reco_pandoraGausCryoE:      [ pandoraGausCryoE,
+icarus_reco_postPandoraGausCryoE:  [ nuslhitsCryoE,
+                                     NuGraphCryoE,
                                      pandoraTrackGausCryoE,
                                      pandoraKalmanTrackGausCryoE,
                                      SBNShowerGausCryoE
@@ -265,6 +279,11 @@ icarus_pandora_Gauss:              [
                                      @sequence::icarus_reco_pandoraGausCryoW
                                    ]
 
+icarus_postPandora_Gauss:          [
+                                     @sequence::icarus_reco_postPandoraGausCryoE,
+                                     @sequence::icarus_reco_postPandoraGausCryoW
+                                   ]
+
 #Add flash matching
 icarus_reco_fm:                    [ fmatchCryoE,
                                      fmatchCryoW,
diff --git a/fcl/reco/Stage1/data/stage1_run2_icarus.fcl b/fcl/reco/Stage1/data/stage1_run2_icarus.fcl
index 53f70ccbd..80219d082 100644
--- a/fcl/reco/Stage1/data/stage1_run2_icarus.fcl
+++ b/fcl/reco/Stage1/data/stage1_run2_icarus.fcl
@@ -6,6 +6,7 @@ physics.reco: [ @sequence::icarus_filter2D_cluster3D,
                 @sequence::icarus_pandora_Gauss,
                 @sequence::icarus_reco_fm,
                 @sequence::icarus_tpcpmtbarycentermatch,
+                @sequence::icarus_postPandora_Gauss,
                 @sequence::icarus_crttrack,
                 @sequence::icarus_crtt0tagging,
                 caloskimCalorimetryCryoE, caloskimCalorimetryCryoW]
diff --git a/fcl/reco/Stage1/mc/stage1_run2_1d_icarus_MC.fcl b/fcl/reco/Stage1/mc/stage1_run2_1d_icarus_MC.fcl
index 85b155397..cbf9f4e40 100644
--- a/fcl/reco/Stage1/mc/stage1_run2_1d_icarus_MC.fcl
+++ b/fcl/reco/Stage1/mc/stage1_run2_1d_icarus_MC.fcl
@@ -6,6 +6,7 @@ physics.reco: [
                 @sequence::icarus_reco_Gauss1D_CryoW ,
                 @sequence::icarus_reco_fm,
                 @sequence::icarus_tpcpmtbarycentermatch,
+                @sequence::icarus_postPandora_Gauss,
                 @sequence::icarus_crttrack,
                 @sequence::icarus_crtt0tagging,
                 caloskimCalorimetryCryoE, caloskimCalorimetryCryoW,
diff --git a/fcl/reco/Stage1/mc/stage1_run2_icarus_MC.fcl b/fcl/reco/Stage1/mc/stage1_run2_icarus_MC.fcl
index 0009c0806..9a1c93ace 100644
--- a/fcl/reco/Stage1/mc/stage1_run2_icarus_MC.fcl
+++ b/fcl/reco/Stage1/mc/stage1_run2_icarus_MC.fcl
@@ -53,6 +53,7 @@ physics.reco: [
                 @sequence::icarus_reco_Gauss2D_CryoW ,
                 @sequence::icarus_reco_fm,
                 @sequence::icarus_tpcpmtbarycentermatch,
+                @sequence::icarus_postPandora_Gauss,
                 @sequence::icarus_crttrack,
                 @sequence::icarus_crtt0tagging,
                 caloskimCalorimetryCryoE, caloskimCalorimetryCryoW,
diff --git a/icaruscode/TPC/CMakeLists.txt b/icaruscode/TPC/CMakeLists.txt
index 3bca0be5a..2652769c2 100644
--- a/icaruscode/TPC/CMakeLists.txt
+++ b/icaruscode/TPC/CMakeLists.txt
@@ -9,3 +9,4 @@ add_subdirectory(Simulation)
 add_subdirectory(Tracking)
 add_subdirectory(ICARUSWireCell)
 add_subdirectory(SPAna)
+add_subdirectory(NuGraph)
diff --git a/icaruscode/TPC/NuGraph/CMakeLists.txt b/icaruscode/TPC/NuGraph/CMakeLists.txt
new file mode 100644
index 000000000..7fdd40479
--- /dev/null
+++ b/icaruscode/TPC/NuGraph/CMakeLists.txt
@@ -0,0 +1,58 @@
+cet_enable_asserts()
+
+set( nugraph_tool_lib_list
+       lardataobj::RecoBase
+       larrecodnn::NuGraphBaseTools
+       TorchScatter::TorchScatter
+       art::Framework_Core
+       hep_concurrency::hep_concurrency
+       art_plugin_types::tool
+)
+
+include_directories($ENV{HEP_HPC_INC})
+
+cet_build_plugin(ICARUSNuGraphLoader art::tool
+  LIBRARIES PRIVATE
+  ${nugraph_tool_lib_list}
+)
+
+simple_plugin(ICARUSNuSliceHitsProducer "module"
+    LIBRARIES PRIVATE
+    art::Framework_Core
+    lardataobj::RecoBase
+    sbnobj::Common_Reco
+)
+
+simple_plugin(ICARUSTrueNuSliceHitsProducer "module"
+    LIBRARIES PRIVATE
+    art::Framework_Core
+    nusimdata::SimulationBase
+    lardataobj::RecoBase
+    lardataalg::DetectorInfo
+    sbnobj::Common_Reco
+    larsim::MCCheater_BackTrackerService_service           
+    larsim::MCCheater_ParticleInventoryService_service
+    lardata::DetectorClocksService
+)
+
+simple_plugin(ICARUSHDF5Maker "module"
+    LIBRARIES PRIVATE
+    art::Framework_Core
+    nusimdata::SimulationBase
+    larsim::MCCheater_BackTrackerService_service           
+    larsim::MCCheater_ParticleInventoryService_service
+    hep_hpc_hdf5
+    lardataobj::RecoBase
+    lardataalg::DetectorInfo
+    lardata::DetectorClocksService
+    larcorealg::Geometry
+    larcore::WireReadout
+    ${HDF5_LIBRARIES}
+)
+
+add_subdirectory(scripts)
+
+install_headers()
+install_fhicl()
+install_source()
+
diff --git a/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc b/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc
new file mode 100644
index 000000000..c6f80bb6d
--- /dev/null
+++ b/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc
@@ -0,0 +1,621 @@
+/**
+ * @file   icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc
+ * @brief  Implementation of `ICARUSHDF5Maker` _art_ module.
+ * @author Giuseppe Cerati (cerati@fnal.gov), V Hewes
+ */
+
+#include "art/Framework/Core/EDAnalyzer.h"
+#include "art/Framework/Core/ModuleMacros.h"
+#include "art/Framework/Principal/Event.h"
+#include "art/Framework/Principal/Handle.h"
+#include "art/Framework/Principal/SubRun.h"
+#include "canvas/Utilities/InputTag.h"
+#include "fhiclcpp/ParameterSet.h"
+#include "canvas/Persistency/Common/FindManyP.h"
+#include "messagefacility/MessageLogger/MessageLogger.h"
+
+#include "larcore/Geometry/Geometry.h"
+#include "larcorealg/Geometry/GeometryCore.h"
+#include "larcorealg/Geometry/PlaneGeo.h"
+#include "larcore/Geometry/WireReadout.h"
+#include "larcorealg/Geometry/WireReadoutGeom.h"
+
+#include "nusimdata/SimulationBase/MCTruth.h"
+#include "lardataobj/RecoBase/Hit.h"
+#include "lardataobj/RecoBase/SpacePoint.h"
+#include "lardataobj/RecoBase/OpHit.h"
+#include "lardataobj/RecoBase/OpFlash.h"
+#include "larcorealg/Geometry/Exceptions.h"
+#include "lardataobj/RecoBase/Slice.h"
+
+#include "larsim/MCCheater/BackTrackerService.h"
+#include "larsim/MCCheater/ParticleInventoryService.h"
+#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
+#include "lardata/DetectorInfoServices/DetectorClocksService.h"
+
+#include "lardataobj/AnalysisBase/BackTrackerMatchingData.h"
+
+#include "hep_hpc/hdf5/make_ntuple.hpp"
+#include "art/Framework/Services/Registry/ServiceHandle.h"
+#include "larcore/CoreUtils/ServiceUtil.h"
+
+#include <set>
+#include <optional>
+
+#include "StitchingUtils.h"
+
+class ICARUSHDF5Maker : public art::EDAnalyzer {
+public:
+  explicit ICARUSHDF5Maker(fhicl::ParameterSet const& p);
+
+  ICARUSHDF5Maker(ICARUSHDF5Maker const&) = delete;
+  ICARUSHDF5Maker(ICARUSHDF5Maker&&) = delete;
+  ICARUSHDF5Maker& operator=(ICARUSHDF5Maker const&) = delete;
+  ICARUSHDF5Maker& operator=(ICARUSHDF5Maker&&) = delete;
+
+  void beginSubRun(art::SubRun const& sr) override;
+  void endSubRun(art::SubRun const& sr) override;
+  void analyze(art::Event const& e) override;
+
+private:
+
+  art::InputTag fTruthLabel;
+  art::InputTag fHitLabel;
+  art::InputTag fHitTruthLabel;
+  art::InputTag fSPLabel;
+  art::InputTag fOpHitLabel;
+  art::InputTag fOpFlashLabel;
+
+  bool fUseMap;
+  std::string fOutputName;
+
+  struct HDFDataFile {
+    hep_hpc::hdf5::File file;  ///< output HDF5 file
+
+    hep_hpc::hdf5::Ntuple<hep_hpc::hdf5::Column<int, 1>,    // event id (run, subrun, event)
+                          hep_hpc::hdf5::Column<int, 1>,    // is cc
+			  hep_hpc::hdf5::Column<int, 1>, // nu pdg
+			  hep_hpc::hdf5::Column<float, 1>,  // nu energy
+			  hep_hpc::hdf5::Column<float, 1>,  // lep energy
+			  hep_hpc::hdf5::Column<float, 1>   // nu dir (x, y, z)
+    > eventNtupleNu; ///< event ntuple with neutrino information
+
+    hep_hpc::hdf5::Ntuple<hep_hpc::hdf5::Column<int, 1>,    // event id (run, subrun, event)
+			  hep_hpc::hdf5::Column<int, 1>,    // spacepoint id
+			  hep_hpc::hdf5::Column<float, 1>,  // 3d position (x, y, z)
+			  hep_hpc::hdf5::Column<int, 1>,     // 2d hit (u, v, y)
+			  hep_hpc::hdf5::Column<float, 1>     //ChiSquared of the hit 
+    > spacePointNtuple; ///< spacepoint ntuple
+
+    hep_hpc::hdf5::Ntuple<hep_hpc::hdf5::Column<int, 1>,    // event id (run, subrun, event)
+    			  hep_hpc::hdf5::Column<int, 1>,    // hit id
+    			  hep_hpc::hdf5::Column<float, 1>,  // hit integral
+    			  hep_hpc::hdf5::Column<float, 1>,  // hit rms
+    			  hep_hpc::hdf5::Column<int, 1>,    // tpc id
+    			  hep_hpc::hdf5::Column<int, 1>,    // global plane
+    			  hep_hpc::hdf5::Column<float, 1>,  // global wire
+    			  hep_hpc::hdf5::Column<float, 1>,  // global time
+    			  hep_hpc::hdf5::Column<int, 1>,    // raw plane
+    			  hep_hpc::hdf5::Column<float, 1>,  // raw wire
+    			  hep_hpc::hdf5::Column<float, 1>,  // raw time
+    			  hep_hpc::hdf5::Column<int, 1>    //cryostat
+    > hitNtuple; ///< hit ntuple
+
+    hep_hpc::hdf5::Ntuple<hep_hpc::hdf5::Column<int, 1>,    // event id (run, subrun, event)
+    			  hep_hpc::hdf5::Column<int, 1>,    // g4 id
+    			  hep_hpc::hdf5::Column<int, 1>,    // particle type
+    			  hep_hpc::hdf5::Column<int, 1>,    // parent g4 id
+    			  hep_hpc::hdf5::Column<int, 1>,    // is from nu
+    			  hep_hpc::hdf5::Column<float, 1>,  // momentum
+    			  hep_hpc::hdf5::Column<float, 1>,  // start position (x, y, z)
+    			  hep_hpc::hdf5::Column<float, 1>,  // end position (x, y, z)
+    			  hep_hpc::hdf5::Column<std::string, 1>, // start process
+    			  hep_hpc::hdf5::Column<std::string, 1>  // end process
+    > particleNtuple; ///< particle ntuple
+
+    hep_hpc::hdf5::Ntuple<hep_hpc::hdf5::Column<int, 1>,    // event id (run, subrun, event)
+    			  hep_hpc::hdf5::Column<int, 1>,    // hit id
+    			  hep_hpc::hdf5::Column<int, 1>,    // g4 id
+    			  hep_hpc::hdf5::Column<float, 1>,  // deposited energy [ MeV ]
+    			  hep_hpc::hdf5::Column<float, 1>,  // x position
+    			  hep_hpc::hdf5::Column<float, 1>,  // y position
+    			  hep_hpc::hdf5::Column<float, 1>   // z position
+    > energyDepNtuple; ///< energy deposition ntuple
+
+    hep_hpc::hdf5::Ntuple<hep_hpc::hdf5::Column<int, 1>,    // event id (run, subrun, event)
+    			  hep_hpc::hdf5::Column<int, 1>,    // hit id
+    			  hep_hpc::hdf5::Column<int, 1>,    // hit_channel
+    			  hep_hpc::hdf5::Column<int, 1>,    // wire pos
+    			  hep_hpc::hdf5::Column<float, 1>,  // peaktime
+    			  hep_hpc::hdf5::Column<float, 1>,  // width
+    			  hep_hpc::hdf5::Column<float, 1>,  // area
+    			  hep_hpc::hdf5::Column<float, 1>,  // amplitude
+    			  hep_hpc::hdf5::Column<float, 1>,  // pe
+    			  hep_hpc::hdf5::Column<int, 1>     // usedInFlash
+    > opHitNtuple; ///< PMT hit ntuple
+
+    hep_hpc::hdf5::Ntuple<hep_hpc::hdf5::Column<int, 1>,    // event id (run, subrun, event)
+    			  hep_hpc::hdf5::Column<int, 1>,    // flash id
+                          hep_hpc::hdf5::Column<int, 1>,    // wire pos
+    			  hep_hpc::hdf5::Column<float, 1>,  // time
+    			  hep_hpc::hdf5::Column<float, 1>,  // time width
+    			  hep_hpc::hdf5::Column<float, 1>,  // Y center
+    			  hep_hpc::hdf5::Column<float, 1>,  // Y width
+    			  hep_hpc::hdf5::Column<float, 1>,  // Z center
+    			  hep_hpc::hdf5::Column<float, 1>,  // Z width
+    			  hep_hpc::hdf5::Column<float, 1>   // totalpe
+    > opFlashNtuple; ///< Flash ntuple
+
+    hep_hpc::hdf5::Ntuple<hep_hpc::hdf5::Column<int, 1>,    // event id (run, subrun, event)
+    			  hep_hpc::hdf5::Column<int, 1>,    // sumpe id
+    			  hep_hpc::hdf5::Column<int, 1>,    // flash id
+    			  hep_hpc::hdf5::Column<int, 1>,    // PMT channel
+                          hep_hpc::hdf5::Column<float, 1>,  // YZ pos
+    			  hep_hpc::hdf5::Column<float, 1>   // pe
+    > opFlashSumPENtuple; ///< Flash SumPE ntuple
+
+    static std::string makeOutputFileName(std::string const& outputName, art::SubRunID const& sr)
+    {
+      struct timeval now;
+      gettimeofday(&now, NULL);
+      std::ostringstream fileName;
+      fileName << outputName
+	       << "_r" << std::setfill('0') << std::setw(5) << sr.run()
+	       << "_s" << std::setfill('0') << std::setw(5) << sr.subRun()
+	       << "_ts" << std::setw(6) << now.tv_usec << ".h5";
+      std::cout << fileName.str() << std::endl;
+      return fileName.str();
+    }
+    
+    HDFDataFile(std::string const& outputName, art::SubRunID const& sr)
+      : file{ makeOutputFileName(outputName, sr), H5F_ACC_TRUNC }
+      , eventNtupleNu{
+      	hep_hpc::hdf5::make_ntuple({file, "event_table", 1000},
+      	hep_hpc::hdf5::make_column<int>("event_id", 3),
+        hep_hpc::hdf5::make_scalar_column<int>("is_cc"),
+        hep_hpc::hdf5::make_scalar_column<int>("nu_pdg"),    
+        hep_hpc::hdf5::make_scalar_column<float>("nu_energy"),
+        hep_hpc::hdf5::make_scalar_column<float>("lep_energy"),
+        hep_hpc::hdf5::make_column<float>("nu_dir", 3))
+        }
+      , spacePointNtuple{
+        hep_hpc::hdf5::make_ntuple({file, "spacepoint_table", 1000},
+        hep_hpc::hdf5::make_column<int>("event_id", 3),
+        hep_hpc::hdf5::make_scalar_column<int>("spacepoint_id"),
+        hep_hpc::hdf5::make_column<float>("position", 3),
+        hep_hpc::hdf5::make_column<int>("hit_id", 3),
+        hep_hpc::hdf5::make_scalar_column<float>("Chi_squared"))
+        }
+      , hitNtuple{
+        hep_hpc::hdf5::make_ntuple({file, "hit_table", 1000},
+        hep_hpc::hdf5::make_column<int>("event_id", 3),
+        hep_hpc::hdf5::make_scalar_column<int>("hit_id"),
+        hep_hpc::hdf5::make_scalar_column<float>("integral"),
+        hep_hpc::hdf5::make_scalar_column<float>("rms"),
+        hep_hpc::hdf5::make_scalar_column<int>("tpc"),
+        hep_hpc::hdf5::make_scalar_column<int>("global_plane"),
+        hep_hpc::hdf5::make_scalar_column<float>("global_wire"),
+        hep_hpc::hdf5::make_scalar_column<float>("global_time"),
+        hep_hpc::hdf5::make_scalar_column<int>("local_plane"),
+        hep_hpc::hdf5::make_scalar_column<float>("local_wire"),
+        hep_hpc::hdf5::make_scalar_column<float>("local_time"),
+      	hep_hpc::hdf5::make_scalar_column<int>("Cryostat"))
+        }
+      , particleNtuple{
+        hep_hpc::hdf5::make_ntuple({file, "particle_table", 1000},
+        hep_hpc::hdf5::make_column<int>("event_id", 3),
+        hep_hpc::hdf5::make_scalar_column<int>("g4_id"),
+        hep_hpc::hdf5::make_scalar_column<int>("type"),
+        hep_hpc::hdf5::make_scalar_column<int>("parent_id"),
+        hep_hpc::hdf5::make_scalar_column<int>("from_nu"),
+        hep_hpc::hdf5::make_scalar_column<float>("momentum"),
+        hep_hpc::hdf5::make_column<float>("start_position", 3),
+        hep_hpc::hdf5::make_column<float>("end_position", 3),
+        hep_hpc::hdf5::make_scalar_column<std::string>("start_process"),
+      	hep_hpc::hdf5::make_scalar_column<std::string>("end_process"))
+        }
+      , energyDepNtuple{
+        hep_hpc::hdf5::make_ntuple({file, "edep_table", 1000},
+        hep_hpc::hdf5::make_column<int>("event_id", 3),
+        hep_hpc::hdf5::make_scalar_column<int>("hit_id"),
+        hep_hpc::hdf5::make_scalar_column<int>("g4_id"),
+        hep_hpc::hdf5::make_scalar_column<float>("energy"),
+        hep_hpc::hdf5::make_scalar_column<float>("x_position"),
+        hep_hpc::hdf5::make_scalar_column<float>("y_position"),
+      	hep_hpc::hdf5::make_scalar_column<float>("z_position"))
+        }
+      , opHitNtuple{
+        hep_hpc::hdf5::make_ntuple({file, "ophit_table", 1000},
+        hep_hpc::hdf5::make_column<int>("event_id", 3),
+        hep_hpc::hdf5::make_scalar_column<int>("hit_id"),
+        hep_hpc::hdf5::make_scalar_column<int>("hit_channel"),
+        hep_hpc::hdf5::make_column<int>("wire_pos", 3),//3 views
+        hep_hpc::hdf5::make_scalar_column<float>("peaktime"),
+        hep_hpc::hdf5::make_scalar_column<float>("width"),
+        hep_hpc::hdf5::make_scalar_column<float>("area"),
+        hep_hpc::hdf5:: make_scalar_column<float>("amplitude"),
+        hep_hpc::hdf5::make_scalar_column<float>("pe"),
+      	hep_hpc::hdf5::make_scalar_column<int>("sumpe_id"))
+        }
+      , opFlashNtuple{
+        hep_hpc::hdf5::make_ntuple({file, "opflash_table", 1000},
+        hep_hpc::hdf5::make_column<int>("event_id", 3),
+        hep_hpc::hdf5::make_scalar_column<int>("flash_id"),
+        hep_hpc::hdf5::make_column<int>("wire_pos", 3),//3 views
+        hep_hpc::hdf5::make_scalar_column<float>("time"),
+        hep_hpc::hdf5::make_scalar_column<float>("time_width"),
+        hep_hpc::hdf5::make_scalar_column<float>("y_center"),
+        hep_hpc::hdf5::make_scalar_column<float>("y_width"),
+        hep_hpc::hdf5::make_scalar_column<float>("z_center"),
+        hep_hpc::hdf5::make_scalar_column<float>("z_width"),
+      	hep_hpc::hdf5::make_scalar_column<float>("totalpe"))
+        }
+      , opFlashSumPENtuple{
+        hep_hpc::hdf5::make_ntuple({file, "opflashsumpe_table", 1000},
+        hep_hpc::hdf5::make_column<int>("event_id", 3),
+        hep_hpc::hdf5::make_scalar_column<int>("sumpe_id"),
+        hep_hpc::hdf5::make_scalar_column<int>("flash_id"),
+        hep_hpc::hdf5::make_scalar_column<int>("pmt_channel"),
+        hep_hpc::hdf5::make_column<float>("yz_pos", 2),
+      	hep_hpc::hdf5::make_scalar_column<float>("sumpe"))
+        }
+      { }
+  };
+  std::unique_ptr<HDFDataFile> fHDFData;
+
+  int NearWire(const geo::PlaneGeo &plane, float x, float y, float z) const;
+};
+
+ICARUSHDF5Maker::ICARUSHDF5Maker(fhicl::ParameterSet const& p)
+  : EDAnalyzer{p},
+    fTruthLabel(p.get<art::InputTag>("TruthLabel")),
+    fHitLabel(  p.get<art::InputTag>("HitLabel")),
+    fHitTruthLabel(  p.get<art::InputTag>("HitTruthLabel","")),
+    fSPLabel(   p.get<art::InputTag>("SPLabel")),
+    fOpHitLabel(  p.get<art::InputTag>("OpHitLabel")),
+    fOpFlashLabel(  p.get<art::InputTag>("OpFlashLabel")),
+    fUseMap(    p.get<bool>("UseMap", false)),
+    fOutputName(p.get<std::string>("OutputName"))
+{ }
+
+void ICARUSHDF5Maker::analyze(art::Event const& e) {
+
+  const cheat::BackTrackerService* bt = fUseMap? nullptr: art::ServiceHandle<cheat::BackTrackerService>().get();
+  geo::WireReadoutGeom const& geom = art::ServiceHandle<geo::WireReadout const>()->Get();
+
+  std::set<int> g4id;
+  auto const* pi = lar::providerFrom<cheat::ParticleInventoryService>();
+  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(e);
+  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(e, clockData);
+  int run = e.id().run();
+  int subrun = e.id().subRun();
+  int event = e.id().event();
+
+  // auto tpcgeo = art::ServiceHandle<geo::Geometry const>()->TPC();
+  // std::cout << "tpc drift distance=" << tpcgeo.DriftDistance() << std::endl;
+  // std::cout << "drift velocity=" << detProp.DriftVelocity() << std::endl;
+  // std::cout << "drift time=" << tpcgeo.DriftDistance()/detProp.DriftVelocity() << " 2*dt/0.4=" << 2*tpcgeo.DriftDistance()/detProp.DriftVelocity()/clockData.TPCClock().TickPeriod() << std::endl;
+  // std::cout << "tick period=" << clockData.TPCClock().TickPeriod() << std::endl;
+  // std::cout << "tpcTime=" << clockData.TPCTime() << std::endl;
+  // std::cout << "triggerOffsetTPC=" << clockData.TriggerOffsetTPC() << std::endl;
+  // std::cout << "triggerTime=" << clockData.TriggerTime() << std::endl;
+  // std::cout << "correction=" << 2*(tpcgeo.DriftDistance()/detProp.DriftVelocity()-clockData.TriggerOffsetTPC())/clockData.TPCClock().TickPeriod() << std::endl;
+
+  // std::cout << "NEW EVENT: " << run <<"  "<< subrun << " "<< event<<std::endl;
+
+  std::array<int, 3> evtID { run, subrun, event };
+  // Fill event table
+  // Get MC truth
+  auto truthHandle = e.getValidHandle<std::vector<simb::MCTruth>>(fTruthLabel);
+  if (truthHandle->size() != 1) {
+    //avoid pile-up, which is not handled downstream
+    return;
+  }
+  simb::MCNeutrino const& nutruth = truthHandle->at(0).GetNeutrino();
+
+  auto up = nutruth.Nu().Momentum().Vect().Unit();
+  std::array<float, 3> nuMomentum {(float)up.X(),(float)up.Y(),(float)up.Z()};
+
+  fHDFData->eventNtupleNu.insert( evtID.data(),
+				  nutruth.CCNC() == simb::kCC,
+				  nutruth.Nu().PdgCode(),
+				  nutruth.Nu().E(),
+				  nutruth.Lepton().E(),
+				  nuMomentum.data()
+				  );
+
+  // for (int ip=0;ip<truthHandle->at(0).NParticles();ip++) {
+  //   std::cout << "mcp tkid=" << truthHandle->at(0).GetParticle(ip).TrackId() << " pdg=" << truthHandle->at(0).GetParticle(ip).PdgCode() 
+  // 		<< " mother=" << truthHandle->at(0).GetParticle(ip).Mother()
+  // 		<< " vtx=" << truthHandle->at(0).GetParticle(ip).Vx() << " " << truthHandle->at(0).GetParticle(ip).Vy() << " " << truthHandle->at(0).GetParticle(ip).Vz()
+  // 		<< std::endl;
+  // }
+
+  mf::LogDebug("ICARUSHDF5Maker") << "Filling event table"
+				  << "\nrun " << evtID[0] << ", subrun " << evtID[1]
+				  << ", event " << evtID[2]
+				  << "\nis cc? " << (nutruth.CCNC() == simb::kCC)
+				  << ", nu energy " << nutruth.Nu().E()
+				  << ", lepton energy " << nutruth.Lepton().E()
+				  << "\nnu momentum x " << nuMomentum[0] << ", y "
+				  << nuMomentum[1] << ", z " << nuMomentum[2];
+
+  // Get spacepoints from the event record
+  auto spListHandle = e.getValidHandle<std::vector<recob::SpacePoint>>(fSPLabel);
+  std::vector<recob::SpacePoint> const& splist = *spListHandle;
+
+  // Get hits from the event record
+  auto hitListHandle = e.getValidHandle<std::vector<recob::Hit>>(fHitLabel);
+  std::vector<art::Ptr<recob::Hit>> hitlist;
+  art::fill_ptr_vector(hitlist, hitListHandle);
+
+  // Get assocations from spacepoints to hits
+  art::FindManyP<recob::Hit> fmp(spListHandle, e, fSPLabel);
+
+  // Fill spacepoint table
+  for (size_t i = 0; i < spListHandle->size(); ++i) {
+
+    std::vector<art::Ptr<recob::Hit>> const& associatedHits = fmp.at(i);
+    if (associatedHits.size() < 3) {
+      throw art::Exception(art::errors::LogicError) << "I am certain this cannot happen... but here you go, space point with " << associatedHits.size() << " hits";
+    }
+
+    std::array<float, 3> pos {(float)splist[i].XYZ()[0],(float)splist[i].XYZ()[1],(float)splist[i].XYZ()[2]};
+
+    std::array<int, 3> hitID { -1, -1, -1 };
+    for (art::Ptr<recob::Hit> const& hit: associatedHits) {
+      int plane = util::stitchedPlane(hit->WireID());
+      //std::cout << " tpc=" << hit->WireID().TPC << " plane=" << hit->WireID().Plane
+      //          << " plane2=" << plane << " key=" <<  hit.key() << std::endl;
+      hitID[plane] = hit.key();
+    }
+    fHDFData->spacePointNtuple.insert(evtID.data(),splist[i].ID(), pos.data(), hitID.data(),splist[i].Chisq() );
+
+    mf::LogDebug("ICARUSHDF5Maker") << "Filling spacepoint table"
+				    << "\nrun " << evtID[0] << ", subrun " << evtID[1]
+				    << ", event " << evtID[2]
+				    << "\nspacepoint id " << splist[i].ID()
+				    << "\nposition x " << pos[0] << ", y " << pos[1]
+				    << ", z " << pos[2]
+				    << "\nhit ids " << hitID[0] << ", " << hitID[1]
+				    << ", " << hitID[2]
+				    << "Chi_squared "  << splist[i].Chisq();
+  } // for spacepoint
+
+  std::unique_ptr<art::FindManyP<simb::MCParticle, anab::BackTrackerHitMatchingData>> hittruth;
+  if (fUseMap) {
+    hittruth = std::unique_ptr<art::FindManyP<simb::MCParticle, anab::BackTrackerHitMatchingData> >(new art::FindManyP<simb::MCParticle, anab::BackTrackerHitMatchingData>(hitListHandle, e, fHitTruthLabel));
+  }
+
+  // Loop over hits
+  for (art::Ptr<recob::Hit> hit : hitlist) {
+
+    // Fill hit table
+    geo::WireID wireid = hit->WireID();
+    size_t plane = util::stitchedPlane(wireid);
+    double time = util::stitchedTime(wireid,hit->PeakTime());
+    size_t wire = util::stitchedWire(wireid);
+    fHDFData->hitNtuple.insert(evtID.data(),
+			       hit.key(), hit->Integral(), hit->RMS(), wireid.TPC,
+			       plane, wire, time,
+			       wireid.Plane, wireid.Wire, hit->PeakTime(),wireid.Cryostat
+			       );
+
+    mf::LogInfo("ICARUSHDF5Maker") << "Filling hit table"
+                             << "\nrun " << evtID[0] << ", subrun " << evtID[1]
+                             << ", event " << evtID[2]
+                             << "\nhit id " << hit.key() << ", integral "
+                             << hit->Integral() << ", RMS " << hit->RMS()
+                             << ", TPC " << wireid.TPC
+                             << "\nglobal plane " << plane << ", global wire "
+                             << wire << ", global time " << time
+                             << "\nlocal plane " << wireid.Plane
+                             << ", local wire " << wireid.Wire
+                             << ", local time " << hit->PeakTime()
+                             << ", Cryostat " << wireid.Cryostat;
+
+    // Fill energy deposit table
+    if (fUseMap) {
+      //not supported in icarus at the moment
+      throw art::Exception{ art::errors::UnimplementedFeature } << "The use of map is currently not supported in ICARUS";
+      /*
+      std::vector<art::Ptr<simb::MCParticle>> particle_vec = hittruth->at(hit.key());
+      std::vector<anab::BackTrackerHitMatchingData const *> match_vec = hittruth->data(hit.key());
+      //loop over particles
+      for (size_t i_p = 0; i_p < particle_vec.size(); ++i_p) {
+        g4id.insert(particle_vec[i_p]->TrackId());
+        fEnergyDepNtuple->insert(evtID.data(), hit.key(),
+                                 particle_vec[i_p]->TrackId(),
+                                 match_vec[i_p]->ideFraction,
+                                 std::numeric_limits<double>::quiet_NaN(),
+                                 std::numeric_limits<double>::quiet_NaN(),
+                                 std::numeric_limits<double>::quiet_NaN());
+        mf::LogInfo("ICARUSHDF5Maker") << "Filling energy deposit table"
+				 << "\nrun " << evtID[0] << ", subrun " << evtID[1]
+				 << ", event " << evtID[2]
+				 << "\nhit id " << hit.key() << ", g4 id "
+				 << particle_vec[i_p]->TrackId() << ", energy fraction "
+				 << match_vec[i_p]->ideFraction << ", position ("
+				 << std::nan << ", " << std::nan << ", " << std::nan << ")";;
+      } // for particle 
+      */
+    }  else {
+
+      std::vector<sim::TrackIDE> const& h_ides = bt->ChannelToTrackIDEs(clockData, hit->Channel(), hit->StartTick(), hit->EndTick());
+      for (auto const& tide : h_ides) {
+      	int tid = tide.trackID;
+	// if negative id, make sure there is a corresponding particle to look for before taking the abs.
+	// This way negative means no associated particle (a convention that can be used in the code that processes the ntuple).
+      	if (pi->TrackIdToParticle_P(abs(tid))) tid = abs(tid);
+        g4id.insert(tid);
+        fHDFData->energyDepNtuple.insert(evtID.data(),hit.key(), tid, tide.energyFrac, -1., -1., -1.);
+      }
+
+    } // if using microboone map method or not
+  } // for hit
+
+  std::unordered_map<int, simb::MCParticle const*> allIDs;
+
+  // Add invisible particles to hierarchy
+  for (int id : g4id) {
+    const simb::MCParticle* p = pi->TrackIdToParticle_P(abs(id));
+    allIDs.emplace(abs(id), p);
+    while (p && p->Mother() != 0 ) {
+      auto mid = abs(p->Mother());
+      p = pi->TrackIdToParticle_P(mid);
+      allIDs.emplace(mid, p);
+    }
+  }
+  // Loop over true particles and fill table
+  for (auto [ id, p ] : allIDs) {
+    if (!p) {
+      mf::LogError("ICARUSHDF5Maker") << "Track ID=" << id << " not tracked back to any particle.";
+      continue;
+    }
+    auto mct = pi->TrackIdToMCTruth_P(abs(id));
+    int fromNu = (mct.isAvailable() ? mct->NeutrinoSet() : 0);
+    // std::cout << "all mcp tkid=" << p->TrackId() << " pdg=" << p->PdgCode()
+    // 	      << " mother=" << p->Mother()
+    // 	      << " vtx=" << p->Vx() << " " << p->Vy() << " " << p->Vz()
+    // 	      << " mctruth=" << mct->NeutrinoSet()
+    // 	      << std::endl;
+    std::array<float, 3> particleStart { (float)p->Vx(), (float)p->Vy(), (float)p->Vz() };
+    std::array<float, 3> particleEnd { (float)p->EndX(), (float)p->EndY(), (float)p->EndZ() };
+    fHDFData->particleNtuple.insert(evtID.data(),
+				    abs(id), p->PdgCode(), p->Mother(), fromNu, (float)p->P(),
+				    particleStart.data(), particleEnd.data(),
+				    p->Process(), p->EndProcess()
+				    );
+    mf::LogDebug("ICARUSHDF5Maker") << "Filling particle table"
+				    << "\nrun " << evtID[0] << ", subrun " << evtID[1]
+				    << ", event " << evtID[2]
+				    << "\ng4 id " << abs(id) << ", pdg code "
+				    << p->PdgCode() << ", parent " << p->Mother()
+				    << ", momentum " << p->P()
+				    << "\nparticle start x " << particleStart[0]
+				    << ", y " << particleStart[1]
+				    << ", z " << particleStart[2]
+				    << "\nparticle end x " << particleEnd[0] << ", y "
+				    << particleEnd[1] << ", z " << particleEnd[2]
+				    << "\nstart process " << p->Process()
+				    << ", end process " << p->EndProcess();
+  }
+
+  // Get OpFlashs from the event record
+  art::Handle< std::vector< recob::OpFlash > > opFlashListHandle;
+  std::optional<art::FindManyP<recob::OpHit> > assocFlashHit;
+  std::vector<recob::OpFlash> const* opflashlist = nullptr;
+  if (!fOpFlashLabel.empty()) {
+    e.getByLabel(fOpFlashLabel, opFlashListHandle);
+    opflashlist = opFlashListHandle.product();
+    assocFlashHit.emplace(opFlashListHandle, e, fOpFlashLabel);
+  }
+
+  // get the flash matching the slice we selected
+  auto const& sliceOpFlashAssns = e.getProduct<art::Assns<recob::Slice, recob::OpFlash>>(fHitLabel);
+  int flkey = sliceOpFlashAssns.size()==0 ? -1: sliceOpFlashAssns.at(0).second.key();
+
+  std::vector<std::vector<int> > flashsumpepmtmap;//this has at most one element, due to the check on flkey
+
+  // Pick only the flash we care about
+  if (flkey>=0) {
+    recob::OpFlash const& opflash = opflashlist->at(flkey);
+    std::vector<float> pes;
+    for (auto pe : opflash.PEs()) pes.push_back(pe);
+    std::vector<int> nearwires;
+    double xyz[3] = {0.,opflash.YCenter(),opflash.ZCenter()};
+    for (auto const& p : geom.Iterate<geo::PlaneGeo>()) nearwires.push_back(NearWire(p,xyz[0],xyz[1],xyz[2]));
+
+    // Fill opflash table
+    fHDFData->opFlashNtuple.insert(evtID.data(),
+				   flkey,
+				   nearwires.data(),
+				   opflash.Time(),opflash.TimeWidth(),
+				   opflash.YCenter(),opflash.YWidth(),opflash.ZCenter(),opflash.ZWidth(),
+				   opflash.TotalPE()
+				   );
+    size_t count = 0;
+    std::vector<int> sumpepmtmap(art::ServiceHandle<geo::Geometry>()->NOpDets(),0);
+    for (size_t ipmt=0;ipmt<pes.size();ipmt++) {
+      if (pes[ipmt]<=0.) continue;
+      auto xyz = geom.OpDetGeoFromOpChannel(ipmt).GetCenter();
+      std::vector<float> yzpos = {float(xyz.Y()),float(xyz.Z())};
+      fHDFData->opFlashSumPENtuple.insert(evtID.data(),count,flkey,ipmt,yzpos.data(),pes[ipmt]);
+      sumpepmtmap[ipmt] = count;
+      count++;
+    }
+    flashsumpepmtmap.push_back(sumpepmtmap);
+    mf::LogDebug("ICARUSHDF5Maker") << "Filling opflash table"
+				    << "\nrun " << evtID[0] << ", subrun " << evtID[1]
+				    << ", event " << evtID[2]
+				    << "\nflash id " << flkey << ", Time " << opflash.Time()
+				    << ", TotalPE " << opflash.TotalPE()//;
+				    << "\nWireCenters size0 " << opflash.WireCenters().size()//;
+				    << "\nYCenter " << opflash.YCenter()<< " ZCenter " << opflash.ZCenter()
+				    << "\nYWidth " << opflash.YWidth()<< " ZWidth " << opflash.ZWidth()
+				    << "\nInBeamFrame " << opflash.InBeamFrame()<< " OnBeamTime " << opflash.OnBeamTime()
+				    << "\nAbsTime " << opflash.AbsTime() << " TimeWidth " << opflash.TimeWidth() << " FastToTotal " << opflash.FastToTotal();
+  }
+
+  // Get OpHits from the event record
+  std::vector< art::Ptr< recob::OpHit > > ophitlist;
+  if (!fOpHitLabel.empty()) {
+    auto opHitListHandle = e.getValidHandle< std::vector< recob::OpHit > >(fOpHitLabel);
+    art::fill_ptr_vector(ophitlist, opHitListHandle);
+  }
+
+  std::set<art::Ptr<recob::OpHit>> ophv;
+  if (flkey >= 0) {
+    auto const& flashHits = assocFlashHit->at(flkey);
+    ophv.insert(begin(flashHits), end(flashHits));
+  }
+  // Loop over ophits
+  for (art::Ptr< recob::OpHit > ophit : ophitlist) {
+    std::vector<int> nearwires;
+    auto xyz = geom.OpDetGeoFromOpChannel(ophit->OpChannel()).GetCenter();
+    for (auto const& p : geom.Iterate<geo::PlaneGeo>()) nearwires.push_back(NearWire(p,xyz.X(),xyz.Y(),xyz.Z()));
+
+    bool isInFlash = (ophv.count(ophit) > 0);
+
+    // Fill ophit table
+    fHDFData->opHitNtuple.insert(evtID.data(),ophit.key(),
+				 ophit->OpChannel(),nearwires.data(),
+				 ophit->PeakTime(),ophit->Width(),
+				 ophit->Area(),ophit->Amplitude(),ophit->PE(),
+				 (isInFlash ? flashsumpepmtmap[0][ophit->OpChannel()] : -1)
+				 );
+    mf::LogDebug("ICARUSHDF5Maker") << "\nFilling ophit table"
+				    << "\nrun " << evtID[0] << ", subrun " << evtID[1]
+				    << ", event " << evtID[2]
+				    << "\nhit id " << ophit.key() << ", channel "
+				    << ophit->OpChannel() << ", PeakTime " << ophit->PeakTime()
+				    << ", Width " << ophit->Width()
+				    << "\narea " << ophit->Area() << ", amplitude "
+				    << ophit->Amplitude() << ", PE " << ophit->PE();
+  }
+  //End optical analyzer
+
+} // function ICARUSHDF5Maker::analyze
+
+void ICARUSHDF5Maker::beginSubRun(art::SubRun const& sr) {
+  fHDFData = std::make_unique<HDFDataFile>(fOutputName, sr.id());
+}
+
+void ICARUSHDF5Maker::endSubRun(art::SubRun const& sr) {
+  fHDFData.reset();
+}
+
+int ICARUSHDF5Maker::NearWire(const geo::PlaneGeo &plane, float x, float y, float z) const
+{
+  geo::WireID wireID;
+  try {
+    wireID = plane.NearestWireID(geo::Point_t(x,y,z));
+  }
+  catch (geo::InvalidWireError const& e) {
+    if (!e.hasSuggestedWire()) throw;
+    wireID = plane.ClosestWireID(e.suggestedWireID());
+  }
+  return wireID.Wire;
+}
+
+DEFINE_ART_MODULE(ICARUSHDF5Maker)
diff --git a/icaruscode/TPC/NuGraph/ICARUSNuGraphLoader_tool.cc b/icaruscode/TPC/NuGraph/ICARUSNuGraphLoader_tool.cc
new file mode 100644
index 000000000..0774c17b5
--- /dev/null
+++ b/icaruscode/TPC/NuGraph/ICARUSNuGraphLoader_tool.cc
@@ -0,0 +1,152 @@
+/**
+ * @file   icaruscode/TPC/NuGraph/ICARUSNuGraphLoader_tool.cc
+ * @author Giuseppe Cerati (cerati@fnal.gov)
+ */
+
+#include "larrecodnn/NuGraph/Tools/LoaderToolBase.h"
+
+#include "art/Utilities/ToolMacros.h"
+
+#include "messagefacility/MessageLogger/MessageLogger.h"
+#include "canvas/Persistency/Common/FindManyP.h"
+#include "canvas/Persistency/Common/Ptr.h"
+#include "canvas/Utilities/InputTag.h"
+#include "lardataobj/RecoBase/SpacePoint.h"
+#include "lardataobj/RecoBase/Hit.h"
+#include "larcoreobj/SimpleTypesAndConstants/geo_types.h" // geo::WireID
+#include <utility> // std::move()
+#include <torch/torch.h>
+
+#include "StitchingUtils.h"
+
+class ICARUSNuGraphLoader : public LoaderToolBase {
+
+/**
+ * @brief Tool to collect the inputs needed by NuGraph from ICARUS data products.
+ *
+ * Read hit and space points from the event record, and package them for usage in NuGraph. This tool is called from larrecodnn/NuGraph/NuGraphInference_module.cc.
+ * Hits are stitched so that the 4 ICARUS TPCs in a cryostat are viewed in a single time vs wire plane.
+ * Only space points with chi2<MinChiSq (currently set to 0.5) and with 3 associated hits are considered.
+ *
+ */
+
+public:
+  /**
+   *  @brief  Constructor
+   *
+   *  @param  pset
+   */
+  ICARUSNuGraphLoader(const fhicl::ParameterSet& pset);
+
+  /**
+   * @brief loadData function
+   *
+   * @param art::Event event record, list of input, idsmap
+   */
+  void loadData(art::Event& e,
+                vector<art::Ptr<recob::Hit>>& hitlist,
+                vector<NuGraphInput>& inputs,
+                vector<vector<size_t>>& idsmap) override;
+
+private:
+
+  art::InputTag const fHitInput;
+  art::InputTag const fSpsInput;
+  static constexpr double MinChiSq = 0.5; ///< Threshold to consider a space point good.
+};
+
+ICARUSNuGraphLoader::ICARUSNuGraphLoader(const fhicl::ParameterSet& p)
+  : fHitInput{p.get<art::InputTag>("hitInput")}, fSpsInput{p.get<art::InputTag>("spsInput")}
+{}
+
+void ICARUSNuGraphLoader::loadData(art::Event& e,
+                              vector<art::Ptr<recob::Hit>>& hitlist,
+                              vector<NuGraphInput>& inputs,
+                              vector<vector<size_t>>& idsmap)
+{
+  //
+  auto hitListHandle = e.getValidHandle<vector<recob::Hit>>(fHitInput);
+  art::fill_ptr_vector(hitlist, hitListHandle);
+  //
+  idsmap = std::vector<std::vector<size_t>>(planes.size(), std::vector<size_t>());
+  for (auto h : hitlist) {
+    //
+    size_t plane = util::stitchedPlane(h->WireID());
+    idsmap[plane].push_back(h.key());
+  }
+
+  vector<int32_t> hit_table_hit_id_data;
+  vector<int32_t> hit_table_local_plane_data;
+  vector<float> hit_table_local_time_data;
+  vector<int32_t> hit_table_local_wire_data;
+  vector<float> hit_table_integral_data;
+  vector<float> hit_table_rms_data;
+  vector<int32_t> spacepoint_table_spacepoint_id_data;
+  vector<int32_t> spacepoint_table_hit_id_u_data;
+  vector<int32_t> spacepoint_table_hit_id_v_data;
+  vector<int32_t> spacepoint_table_hit_id_y_data;
+
+  // hit table
+  for (auto h : hitlist) {
+
+    geo::WireID wireid = h->WireID();
+    size_t plane = util::stitchedPlane(wireid);
+    double time = util::stitchedTime(wireid,h->PeakTime());
+    size_t wire = util::stitchedWire(wireid);
+
+    hit_table_hit_id_data.push_back(h.key());
+    hit_table_local_plane_data.push_back(plane);
+    hit_table_local_time_data.push_back(time);
+    hit_table_local_wire_data.push_back(wire);
+    hit_table_integral_data.push_back(h->Integral());
+    hit_table_rms_data.push_back(h->RMS());
+  }
+  mf::LogDebug{ "ICARUSNuGraphLoader" } << "loader has nhits=" << hit_table_hit_id_data.size();
+
+  // Get spacepoints from the event record
+  art::Handle<vector<recob::SpacePoint>> spListHandle;
+  vector<art::Ptr<recob::SpacePoint>> splist;
+  if (e.getByLabel(fSpsInput, spListHandle)) { art::fill_ptr_vector(splist, spListHandle); }
+  // Get assocations from spacepoints to hits
+  vector<vector<art::Ptr<recob::Hit>>> sp2Hit(splist.size());
+  if (!splist.empty()) {
+    art::FindManyP<recob::Hit> fmp(spListHandle, e, fSpsInput);
+    for (size_t spIdx = 0; spIdx < sp2Hit.size(); ++spIdx) {
+      if (splist[spIdx]->Chisq()>MinChiSq) continue;
+      // only space points with hits on all planes are enough for NuGraph
+      if (fmp.at(spIdx).size()!=3) continue;
+      sp2Hit[spIdx] = fmp.at(spIdx);
+    }
+  }
+
+  // space point table
+  size_t spidx = 0;
+  for (size_t i = 0; i < splist.size(); ++i) {
+    if (sp2Hit[i].empty()) continue;
+    spacepoint_table_spacepoint_id_data.push_back(spidx++);
+    spacepoint_table_hit_id_u_data.push_back(-1);
+    spacepoint_table_hit_id_v_data.push_back(-1);
+    spacepoint_table_hit_id_y_data.push_back(-1);
+    for (size_t j = 0; j < sp2Hit[i].size(); ++j) {
+      //
+      size_t plane = util::stitchedPlane(sp2Hit[i][j]->WireID());
+      if (plane == 0) spacepoint_table_hit_id_u_data.back() = sp2Hit[i][j].key();
+      if (plane == 1) spacepoint_table_hit_id_v_data.back() = sp2Hit[i][j].key();
+      if (plane == 2) spacepoint_table_hit_id_y_data.back() = sp2Hit[i][j].key();
+    }
+  }
+  mf::LogDebug{ "ICARUSNuGraphLoader" }  << "loader has nsps=" << spacepoint_table_hit_id_u_data.size();
+
+  inputs.emplace_back("hit_table_hit_id", std::move(hit_table_hit_id_data));
+  inputs.emplace_back("hit_table_local_plane", std::move(hit_table_local_plane_data));
+  inputs.emplace_back("hit_table_local_time", std::move(hit_table_local_time_data));
+  inputs.emplace_back("hit_table_local_wire", std::move(hit_table_local_wire_data));
+  inputs.emplace_back("hit_table_integral", std::move(hit_table_integral_data));
+  inputs.emplace_back("hit_table_rms", std::move(hit_table_rms_data));
+
+  inputs.emplace_back("spacepoint_table_spacepoint_id", std::move(spacepoint_table_spacepoint_id_data));
+  inputs.emplace_back("spacepoint_table_hit_id_u", std::move(spacepoint_table_hit_id_u_data));
+  inputs.emplace_back("spacepoint_table_hit_id_v", std::move(spacepoint_table_hit_id_v_data));
+  inputs.emplace_back("spacepoint_table_hit_id_y", std::move(spacepoint_table_hit_id_y_data));
+}
+DEFINE_ART_CLASS_TOOL(ICARUSNuGraphLoader)
diff --git a/icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc b/icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc
new file mode 100644
index 000000000..3b48cb6d6
--- /dev/null
+++ b/icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc
@@ -0,0 +1,187 @@
+/**
+ * @file   icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc
+ * @brief  Implementation of `ICARUSNuSliceHitsProducer` _art_ module.
+ * @author Giuseppe Cerati (cerati@fnal.gov)
+ * @date   May 25, 2021
+ */
+
+#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 "canvas/Utilities/InputTag.h"
+#include "fhiclcpp/ParameterSet.h"
+#include "messagefacility/MessageLogger/MessageLogger.h"
+#include "canvas/Persistency/Common/FindOneP.h"
+#include "canvas/Persistency/Common/FindManyP.h"
+
+#include <algorithm> // std::find()
+#include <cmath> // std::abs()
+#include <memory>
+#include <unordered_map>
+#include <utility> // std::move()
+#include <vector>
+
+#include "lardataobj/RecoBase/OpFlash.h"
+#include "canvas/Persistency/Common/Assns.h"
+#include "art/Persistency/Common/PtrMaker.h"
+#include "lardataobj/RecoBase/Hit.h"
+#include "lardataobj/RecoBase/Slice.h"
+#include "lardataobj/RecoBase/SpacePoint.h"
+
+#include "sbnobj/Common/Reco/TPCPMTBarycenterMatch.h"
+
+class ICARUSNuSliceHitsProducer : public art::EDProducer {
+
+/**
+ * @brief Produce separate collections of hits and space point from the slice identified as most likely from the neutrino interaction.
+ *
+ * Produce separate collections of hits and space point from the slice identified as most likely from the neutrino interaction.
+ * The slice is the one identified as the one producing the trigger flash. Only space points made out of 3 hits are saved.
+ * It also produces a vector of ints, that maps the new hit collection to the original one.
+ * Optionally, it can store the association to the trigger flash used to pick the slice.
+ *
+ */
+
+public:
+  explicit ICARUSNuSliceHitsProducer(fhicl::ParameterSet const& p);
+  // The compiler-generated destructor is fine for non-base
+  // classes without bare pointers or other resource use.
+
+  // Plugins should not be copied or assigned.
+  ICARUSNuSliceHitsProducer(ICARUSNuSliceHitsProducer const&) = delete;
+  ICARUSNuSliceHitsProducer(ICARUSNuSliceHitsProducer&&) = delete;
+  ICARUSNuSliceHitsProducer& operator=(ICARUSNuSliceHitsProducer const&) = delete;
+  ICARUSNuSliceHitsProducer& operator=(ICARUSNuSliceHitsProducer&&) = delete;
+
+
+private:
+
+  // Declare member data here.
+  art::InputTag const fBaryMatchLabel;
+  art::InputTag const fSliceLabel;
+  art::InputTag const fSpsLabel;
+  bool const doFlashAssn;
+
+  // Required functions.
+  void produce(art::Event& e) override;
+
+};
+
+
+ICARUSNuSliceHitsProducer::ICARUSNuSliceHitsProducer(fhicl::ParameterSet const& p)
+  : EDProducer{p},
+  fBaryMatchLabel( p.get<art::InputTag>("BaryMatchLabel","tpcpmtbarycentermatchCryoE")),
+  fSliceLabel(p.get<art::InputTag>("SliceLabel","pandoraGausCryoE")),
+  fSpsLabel(p.get<art::InputTag>("SpsLabel","cluster3DCryoE")),
+  doFlashAssn(p.get<bool>("DoFlashAssn",false))
+  // More initializers here.
+{
+  // Call appropriate produces<>() functions here.
+  produces<std::vector<recob::Hit>>();
+  produces<std::vector<unsigned int>>();
+  produces<std::vector<recob::SpacePoint>>();
+  produces<art::Assns<recob::Hit, recob::SpacePoint>>();
+  if (doFlashAssn) produces<art::Assns<recob::Slice, recob::OpFlash>>();
+
+  // Call appropriate consumes<>() for any products to be retrieved by this module.
+}
+
+void ICARUSNuSliceHitsProducer::produce(art::Event& e)
+{
+
+  auto outputHits = std::make_unique<std::vector<recob::Hit> >();
+  auto assocSliceHitKeys = std::make_unique<std::vector<unsigned int> >();
+
+  art::PtrMaker<recob::Hit> hitPtrMaker(e);
+  art::PtrMaker<recob::SpacePoint> spsPtrMaker(e);
+
+  auto outputSpacepoints = std::make_unique<std::vector<recob::SpacePoint> >();
+  auto outputSpsHitAssns = std::make_unique<art::Assns<recob::Hit, recob::SpacePoint>>();
+  auto outputSlcFlhAssns = std::make_unique<art::Assns<recob::Slice, recob::OpFlash>>();
+
+  art::ValidHandle< std::vector< sbn::TPCPMTBarycenterMatch > > baryMatchListHandle = e.getValidHandle<std::vector<sbn::TPCPMTBarycenterMatch> >(fBaryMatchLabel);
+  art::FindOneP<recob::Slice> msl(baryMatchListHandle, e, fBaryMatchLabel);
+  art::FindOneP<recob::OpFlash> mfl(baryMatchListHandle, e, fBaryMatchLabel);
+  float minRad = 99999.;
+  art::Ptr<recob::Slice> triggeringSlice;
+  art::Ptr<recob::OpFlash> triggeringFlash;
+  for (size_t ibm=0; ibm<baryMatchListHandle->size(); ++ibm) {
+    sbn::TPCPMTBarycenterMatch const& bm = (*baryMatchListHandle)[ibm];
+    // require the radius to be set (i.e. that there is a best match) and to match the one of the trigger within an arbitrarily small enough margin
+    if (bm.radius_Trigger<=0 || std::abs(bm.radius-bm.radius_Trigger)>=0.00001) continue;
+    art::Ptr<recob::Slice> const& slice = msl.at(ibm);
+    art::Ptr<recob::OpFlash> const& flash = mfl.at(ibm);
+    if (!slice || !flash) throw std::logic_error("Unexpected missing flash or slice in barycenter matching associations");
+    mf::LogTrace{ "ICARUSNuSliceHitsProducer" }
+      << "ibm=" << ibm << " radius=" << bm.radius << " radius_Trigger=" << bm.radius_Trigger << " flashTime=" << bm.flashTime
+      << " slkey=" << msl.at(ibm).key() << " center=" << bm.chargeCenter
+      << " flkey=" << mfl.at(ibm).key() << " center=" << bm.flashCenter;
+    if (bm.radius_Trigger<minRad) {
+      minRad = bm.radius_Trigger;
+      triggeringSlice = slice;
+      triggeringFlash = flash;
+    }
+  }
+  mf::LogTrace{ "ICARUSNuSliceHitsProducer" } << "keys: slice=" << triggeringSlice.key() << " flash=" << triggeringFlash.key();
+
+  if (doFlashAssn) {
+    if (triggeringSlice && triggeringFlash) {
+      outputSlcFlhAssns->addSingle(triggeringSlice,triggeringFlash);
+    }
+  }
+
+  art::ValidHandle<std::vector<recob::Slice> > inputSlice = e.getValidHandle<std::vector<recob::Slice> >(fSliceLabel);
+  art::FindManyP<recob::Hit> assocSliceHit(inputSlice, e, fSliceLabel);
+
+  if (triggeringSlice) {
+    mf::LogTrace{ "ICARUSNuSliceHitsProducer" } << "slice hits=" << assocSliceHit.at(triggeringSlice.key()).size();
+    auto const& hits = assocSliceHit.at(triggeringSlice.key());
+    std::unordered_map<size_t,size_t> keyIdxMap;
+    for (size_t ih=0; ih<hits.size(); ih++) {
+      assocSliceHitKeys->push_back(hits[ih].key());
+      keyIdxMap.insert({hits[ih].key(), ih});
+    }
+    std::sort(begin(*assocSliceHitKeys), end(*assocSliceHitKeys));
+    outputHits->reserve(assocSliceHitKeys->size());
+    for (auto const& key: *assocSliceHitKeys) outputHits->push_back(*hits[keyIdxMap.at(key)]);
+  }
+
+  // Get spacepoints from the event record
+  auto splist = e.getValidHandle<std::vector<recob::SpacePoint>>(fSpsLabel);
+  size_t countsps = 0;
+  // Get assocations from spacepoints to hits
+  art::FindManyP<recob::Hit> spacePointToHits(splist, e, fSpsLabel);
+  for (size_t spIdx = 0; spIdx < splist->size(); ++spIdx) {
+    auto const& assochits = spacePointToHits.at(spIdx);
+    // Consider only space points with hits associated on all planes. This is enough for NuGraph.
+    if (assochits.size()<3) continue;
+    //if (splist[spIdx]->Chisq()>0.5) continue;
+    //
+    std::vector<size_t> hitpos;
+    for (size_t j = 0; j < assochits.size(); ++j) {
+      auto pos = std::lower_bound(assocSliceHitKeys->begin(),assocSliceHitKeys->end(),assochits[j].key());
+      if ( (pos == assocSliceHitKeys->end()) || (*pos != assochits[j].key()) ) break;
+      hitpos.push_back(pos-assocSliceHitKeys->begin());
+    }
+    if (hitpos.size()<3) continue;
+    outputSpacepoints->emplace_back((*splist)[spIdx]);
+    //
+    const art::Ptr<recob::SpacePoint> sps = spsPtrMaker(outputSpacepoints->size()-1);
+    for (size_t j = 0; j < hitpos.size(); ++j) {
+      const art::Ptr<recob::Hit> ahp = hitPtrMaker(hitpos[j]);
+      outputSpsHitAssns->addSingle(ahp,sps);
+    }
+    countsps++;
+  } // for spacepoint
+  mf::LogTrace{ "ICARUSNuSliceHitsProducer" } << "sps count=" << countsps;
+
+  e.put(std::move(outputHits));
+  e.put(std::move(assocSliceHitKeys));
+  e.put(std::move(outputSpacepoints));
+  e.put(std::move(outputSpsHitAssns));
+  if (doFlashAssn) e.put(std::move(outputSlcFlhAssns));
+
+}
+
+DEFINE_ART_MODULE(ICARUSNuSliceHitsProducer)
diff --git a/icaruscode/TPC/NuGraph/ICARUSTrueNuSliceHitsProducer_module.cc b/icaruscode/TPC/NuGraph/ICARUSTrueNuSliceHitsProducer_module.cc
new file mode 100644
index 000000000..1b85b3db1
--- /dev/null
+++ b/icaruscode/TPC/NuGraph/ICARUSTrueNuSliceHitsProducer_module.cc
@@ -0,0 +1,210 @@
+/**
+ * @file   icaruscode/TPC/NuGraph/ICARUSTrueNuSliceHitsProducer_module.cc
+ * @brief  Implementation of `ICARUSTrueNuSliceHitsProducer` module
+ * @author Giuseppe Cerati (cerati@fnal.gov)
+ */
+
+#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 "canvas/Utilities/InputTag.h"
+#include "fhiclcpp/ParameterSet.h"
+#include "messagefacility/MessageLogger/MessageLogger.h"
+#include "canvas/Persistency/Common/FindManyP.h"
+
+#include <algorithm> // std::find()
+#include <cmath> // std::abs()
+#include <memory>
+#include <unordered_map>
+#include <utility> // std::move()
+#include <vector>
+
+#include "lardataobj/RecoBase/OpFlash.h"
+#include "canvas/Persistency/Common/Assns.h"
+#include "art/Persistency/Common/PtrMaker.h"
+#include "lardataobj/RecoBase/Hit.h"
+#include "lardataobj/RecoBase/Slice.h"
+#include "lardataobj/RecoBase/PFParticle.h"
+#include "lardataobj/RecoBase/SpacePoint.h"
+
+#include "sbnobj/Common/Reco/TPCPMTBarycenterMatch.h"
+#include "larsim/MCCheater/ParticleInventoryService.h"
+#include "larsim/MCCheater/BackTrackerService.h"
+#include "lardata/DetectorInfoServices/DetectorClocksService.h"
+#include "larcore/CoreUtils/ServiceUtil.h"
+
+class ICARUSTrueNuSliceHitsProducer : public art::EDProducer {
+
+/**
+ * @brief Produce separate collections of hits and space point from the slice with most true hits from the neutrino interaction.
+ *
+ * Produce separate collections of hits and space point from the slice with most true hits from the neutrino interaction.
+ * It also produces a vector of ints, that maps the new hit collection to the original one.
+ * Optionally, it can store the association to the trigger flash matched to the slice.
+ *
+ */
+
+public:
+  explicit ICARUSTrueNuSliceHitsProducer(fhicl::ParameterSet const& p);
+  // The compiler-generated destructor is fine for non-base
+  // classes without bare pointers or other resource use.
+
+  // Plugins should not be copied or assigned.
+  ICARUSTrueNuSliceHitsProducer(ICARUSTrueNuSliceHitsProducer const&) = delete;
+  ICARUSTrueNuSliceHitsProducer(ICARUSTrueNuSliceHitsProducer&&) = delete;
+  ICARUSTrueNuSliceHitsProducer& operator=(ICARUSTrueNuSliceHitsProducer const&) = delete;
+  ICARUSTrueNuSliceHitsProducer& operator=(ICARUSTrueNuSliceHitsProducer&&) = delete;
+
+
+private:
+
+  // Declare member data here.
+  art::InputTag const fBaryMatchLabel;
+  art::InputTag const fSliceLabel;
+  art::InputTag const fSpsLabel;
+  bool const doFlashAssn;
+
+  // Required functions.
+  void produce(art::Event& e) override;
+
+};
+
+
+ICARUSTrueNuSliceHitsProducer::ICARUSTrueNuSliceHitsProducer(fhicl::ParameterSet const& p)
+  : EDProducer{p},
+  fBaryMatchLabel( p.get<art::InputTag>("BaryMatchLabel") ),
+  fSliceLabel(p.get<art::InputTag>("SliceLabel")),
+  fSpsLabel(p.get<art::InputTag>("SpsLabel")),
+  doFlashAssn(p.get<bool>("DoFlashAssn",false))
+  // More initializers here.
+{
+  // Call appropriate produces<>() functions here.
+  produces<std::vector<recob::Hit>>();
+  produces<std::vector<unsigned int>>();
+  produces<std::vector<recob::SpacePoint>>();
+  produces<art::Assns<recob::Hit, recob::SpacePoint>>();
+  if (doFlashAssn) produces<art::Assns<recob::Slice, recob::OpFlash>>();
+
+  // Call appropriate consumes<>() for any products to be retrieved by this module.
+}
+
+void ICARUSTrueNuSliceHitsProducer::produce(art::Event& e)
+{
+
+  auto outputHits = std::make_unique<std::vector<recob::Hit> >();
+  auto assocSliceHitKeys = std::make_unique<std::vector<unsigned int> >();
+
+  art::PtrMaker<recob::Hit> hitPtrMaker(e);
+  art::PtrMaker<recob::SpacePoint> spsPtrMaker(e);
+
+  auto outputSpacepoints = std::make_unique<std::vector<recob::SpacePoint> >();
+  auto outputSpsHitAssns = std::make_unique<art::Assns<recob::Hit, recob::SpacePoint>>();
+  auto outputSlcFlhAssns = std::make_unique<art::Assns<recob::Slice, recob::OpFlash>>();
+
+  auto const* pi = lar::providerFrom<cheat::ParticleInventoryService>();
+  auto const* bt = lar::providerFrom<cheat::BackTrackerService>();
+  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(e);
+
+  art::ValidHandle<std::vector<recob::Slice> > inputSlice = e.getValidHandle<std::vector<recob::Slice> >(fSliceLabel);
+  art::FindManyP<recob::Hit> assocSliceHit(inputSlice, e, fSliceLabel);
+  int slkey = -1;
+  int maxCnt = 0;
+  // find the slice with the largest number of hits contributed by neutrino interactions
+  for (size_t sk=0; sk<inputSlice->size(); sk++) {
+    int slhcnt = 0;
+    for (auto hit : assocSliceHit.at(sk)) {
+      std::vector<sim::TrackIDE> const& h_ides = bt->ChannelToTrackIDEs(clockData, hit->Channel(), hit->StartTick(), hit->EndTick());
+      for (auto const& tide : h_ides) {
+      	int tid = tide.trackID;
+	auto mct = pi->TrackIdToMCTruth_P(abs(tid));
+	bool fromNu = mct.isAvailable() && mct->NeutrinoSet();
+	if (fromNu) {
+	  slhcnt++;
+	  break;
+	}
+      }
+    }
+    if (slhcnt>maxCnt) {
+      slkey = sk;
+      maxCnt = slhcnt;
+    }
+  }
+  mf::LogTrace{ "ICARUSTrueNuSliceHitsProducer" } << "keys: slice=" << slkey;
+
+  if (slkey>=0) {
+    auto const& hits = assocSliceHit.at(slkey);
+    mf::LogTrace{ "ICARUSTrueNuSliceHitsProducer" } << "slice hits=" << hits.size();
+    std::unordered_map<size_t,size_t> keyIdxMap;
+    for (size_t ih=0; ih<hits.size(); ih++) {
+      assocSliceHitKeys->push_back(hits[ih].key());
+      keyIdxMap.insert({hits[ih].key(), ih});
+    }
+    std::sort(begin(*assocSliceHitKeys), end(*assocSliceHitKeys));
+    outputHits->reserve(assocSliceHitKeys->size());
+    for (auto const& key: *assocSliceHitKeys) outputHits->push_back(*hits[keyIdxMap.at(key)]);
+  }
+
+  // Get spacepoints from the event record
+  auto splist = e.getValidHandle<std::vector<recob::SpacePoint>>(fSpsLabel);
+  size_t countsps = 0;
+  // Get assocations from spacepoints to hits
+  art::FindManyP<recob::Hit> spacePointToHits(splist, e, fSpsLabel);
+  for (size_t spIdx = 0; spIdx < splist->size(); ++spIdx) {
+    auto const& assochits = spacePointToHits.at(spIdx);
+    // Consider only space points with hits associated on all planes. This is enough for NuGraph.
+    if (assochits.size()<3) continue;
+    //
+    std::vector<size_t> hitpos;
+    for (size_t j = 0; j < assochits.size(); ++j) {
+      auto pos = std::lower_bound(assocSliceHitKeys->begin(),assocSliceHitKeys->end(),assochits[j].key());
+      if ( (pos == assocSliceHitKeys->end()) || (*pos != assochits[j].key()) ) break;
+      hitpos.push_back(pos-assocSliceHitKeys->begin());
+    }
+    if (hitpos.size()<3) continue;
+    outputSpacepoints->emplace_back((*splist)[spIdx]);
+    //
+    const art::Ptr<recob::SpacePoint> sps = spsPtrMaker(outputSpacepoints->size()-1);
+    for (size_t j = 0; j < hitpos.size(); ++j) {
+      const art::Ptr<recob::Hit> ahp = hitPtrMaker(hitpos[j]);
+      outputSpsHitAssns->addSingle(ahp,sps);
+    }
+    countsps++;
+  } // for spacepoint
+  mf::LogTrace{ "ICARUSTrueNuSliceHitsProducer" } << "sps count=" << countsps;
+
+  art::ValidHandle< std::vector< sbn::TPCPMTBarycenterMatch > > baryMatchListHandle = e.getValidHandle<std::vector<sbn::TPCPMTBarycenterMatch> >(fBaryMatchLabel);
+  art::FindManyP<recob::Slice> msl(baryMatchListHandle, e, fBaryMatchLabel);
+  art::FindManyP<recob::OpFlash> mfl(baryMatchListHandle, e, fBaryMatchLabel);
+  float minRad = 99999.;
+  int mbkey = -1;
+  for (size_t ibm=0; ibm<baryMatchListHandle->size(); ++ibm) {
+    art::Ptr<sbn::TPCPMTBarycenterMatch> bm(baryMatchListHandle,ibm);
+    if (int(msl.at(ibm).at(0).key())!=slkey) continue;
+    if (mfl.at(ibm).size()>0) {
+      mf::LogTrace{ "ICARUSTrueNuSliceHitsProducer" } << "ibm=" << ibm << " radius=" << bm->radius << " radius_Trigger=" << bm->radius_Trigger << " flashTime=" << bm->flashTime
+						  << " slkey=" << msl.at(ibm).at(0).key() << " center=" << bm->chargeCenter
+						  << " flkey=" << mfl.at(ibm).at(0).key() << " center=" << bm->flashCenter;
+      if (bm->radius<minRad) {
+	//this is the triggering flash
+	mbkey = ibm;
+	minRad = bm->radius;
+      }
+    }
+  }
+  mf::LogTrace{ "ICARUSTrueNuSliceHitsProducer" } << "mbkey=" << mbkey;
+  if (doFlashAssn){
+    if (slkey>=0 && mbkey>=0) {
+      art::Ptr<recob::Slice> slp(inputSlice,slkey);
+      outputSlcFlhAssns->addSingle(slp,mfl.at(mbkey).at(0));
+    }
+  }
+
+  e.put(std::move(outputHits));
+  e.put(std::move(assocSliceHitKeys));
+  e.put(std::move(outputSpacepoints));
+  e.put(std::move(outputSpsHitAssns));
+  if (doFlashAssn) e.put(std::move(outputSlcFlhAssns));
+}
+
+DEFINE_ART_MODULE(ICARUSTrueNuSliceHitsProducer)
diff --git a/icaruscode/TPC/NuGraph/StitchingUtils.h b/icaruscode/TPC/NuGraph/StitchingUtils.h
new file mode 100644
index 000000000..b358b6331
--- /dev/null
+++ b/icaruscode/TPC/NuGraph/StitchingUtils.h
@@ -0,0 +1,50 @@
+/**
+ * @file   icaruscode/NuGraph/StitchingUtils.h
+ * @brief  Utilities to stitch together the different logical TPCs in each cryostat.
+ * @author Giuseppe Cerati (cerati@fnal.gov)
+ * @date   June 05, 2025
+ *
+ * This library provides utilities to stitch together the different logical TPCs in each cryostat.
+ * As a result, a single time vs wire coordinate system can be used in each cryostat.
+ *
+ * This is a header-only library.
+ */
+
+#ifndef ICARUSCODE_NUGRAPH_STITCHINGUTILS_H
+#define ICARUSCODE_NUGRAPH_STITCHINGUTILS_H
+
+#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
+
+namespace util {
+
+  int stitchedPlane(const geo::WireID wid) {
+    int plane = wid.Plane;
+    if(wid.TPC==2 || wid.TPC==3) {
+      if(wid.Plane==1) plane=2;
+      else if(wid.Plane==2) plane=1;
+    }
+    return plane;
+  }
+
+  float stitchedTime(const geo::WireID wid, float timein) {
+    float time = timein;
+    if(wid.TPC==2 || wid.TPC==3) {
+      //correction = 2*(tpcgeo.DriftDistance()/detProp.DriftVelocity()-clockData.TriggerOffsetTPC())/clockData.TPCClock().TickPeriod() = 6442.15 us
+      time = 6442.15 - timein;
+    }
+    return time;
+  }
+
+  size_t stitchedWire(const geo::WireID wid) {
+    size_t wire = wid.Wire;
+    int plane = stitchedPlane(wid);
+    if(wid.TPC==1 || wid.TPC==3) {
+      if(plane==1 || plane == 2) {
+	wire = wid.Wire + 2535; //2535 is the last part of the wires in cryos 0 an 2 before the cut in z=0
+      }
+    }
+    return wire;
+  }
+
+} // namespace util
+#endif
diff --git a/icaruscode/TPC/NuGraph/hdf5maker_icarus_slice.fcl b/icaruscode/TPC/NuGraph/hdf5maker_icarus_slice.fcl
new file mode 100644
index 000000000..399d72c6a
--- /dev/null
+++ b/icaruscode/TPC/NuGraph/hdf5maker_icarus_slice.fcl
@@ -0,0 +1,94 @@
+###
+### fcl file for producing the hdf5 files used as input for NuGraph training.
+### Author: G. Cerati (FNAL)
+###
+
+#include "services_common_icarus.fcl"
+#include "services_icarus_simulation.fcl"
+
+process_name: makehdf5
+
+services:
+{
+  TFileService:          { fileName: "reco_hist.root" }
+  TimeTracker:           {}
+  MemoryTracker:         {}
+  RandomNumberGenerator: {}
+  FileCatalogMetadata:   @local::art_file_catalog_mc
+  @table::icarus_common_services
+  @table::icarus_backtracking_services 
+  message:               @local::icarus_message_services_prod_debug
+  ParticleInventoryService: { 
+        ParticleInventory: {
+        EveIdCalculator: "EmEveIdCalculator"
+        G4ModuleLabel: "largeant"
+	OverrideRealData: true
+     }
+  }
+
+                         
+} # services
+
+source:
+{
+  module_type: RootInput
+  maxEvents:   -1
+}
+
+physics:
+{
+
+  analyzers:
+  {
+    hdf5maker: {
+      module_type:  "ICARUSHDF5Maker"
+      TruthLabel:   "generator"
+      SPLabel: "nuslhits"
+      HitLabel: "nuslhits"
+      OpFlashLabel:"opflashCryoE"
+      OpHitLabel:"ophit"
+      UseMap: false
+      OutputName:   "NeutrinoML"
+    }
+  }
+
+ producers:
+  {
+    #nuslhitsE: { 
+    #  module_type: "ICARUSNuSliceHitsProducer"
+    #  BaryMatchLabel: "tpcpmtbarycentermatchCryoE"
+    #  SliceLabel:    "pandoraGausCryoE"
+    #  SpsLabel:      "cluster3DCryoE"
+    #}
+    #nuslhitsW: { 
+    #  module_type: "ICARUSNuSliceHitsProducer"
+    #  BaryMatchLabel: "tpcpmtbarycentermatchCryoW"
+    #  SliceLabel:    "pandoraGausCryoW"
+    #  SpsLabel:      "cluster3DCryoW"
+    #}
+    #nuslhits: {
+    #  module_type:             "HitMerger"
+    #  HitProducerLabelVec:     ["nuslhitsE","nuslhitsW"]
+    #}
+    nuslhits: { 
+      #module_type: "ICARUSNuSliceHitsProducer"
+      module_type: "ICARUSTrueNuSliceHitsProducer"
+      BaryMatchLabel: "tpcpmtbarycentermatchCryoE"
+      SliceLabel:    "pandoraGausCryoE"
+      SpsLabel:      "cluster3DCryoE"
+      DoFlashAssn: true
+    }
+    # sps: cluster3DCryoE
+  }
+ 
+ #reco: [nuslhitsE,nuslhitsW,nuslhits]
+ reco: [nuslhits]
+ ana:       [ hdf5maker ]
+  
+ trigger_paths:       [ reco ]
+ end_paths:     [ ana ]
+  
+}
+
+services.BackTrackerService.BackTracker.G4ModuleLabel: "largeant"
+services.BackTrackerService.BackTracker.SimChannelModuleLabel: "merge"
diff --git a/icaruscode/TPC/NuGraph/nugraph_icarus.fcl b/icaruscode/TPC/NuGraph/nugraph_icarus.fcl
new file mode 100644
index 000000000..d9adcd3aa
--- /dev/null
+++ b/icaruscode/TPC/NuGraph/nugraph_icarus.fcl
@@ -0,0 +1,38 @@
+#include "nugraph.fcl"
+
+BEGIN_PROLOG
+
+nuslhitsCryoE: { 
+  module_type:    "ICARUSNuSliceHitsProducer"
+  BaryMatchLabel: "tpcpmtbarycentermatchCryoE"
+  SliceLabel:     "pandoraGausCryoE"
+  SpsLabel:       "cluster3DCryoE"
+}
+
+nuslhitsCryoW: { 
+  module_type:    "ICARUSNuSliceHitsProducer"
+  BaryMatchLabel: "tpcpmtbarycentermatchCryoW"
+  SliceLabel:     "pandoraGausCryoW"
+  SpsLabel:       "cluster3DCryoW"
+}
+
+NuGraphCryoE: @local::NuGraphLibTorch
+NuGraphCryoE.LoaderTool.tool_type: "ICARUSNuGraphLoader"
+NuGraphCryoE.LoaderTool.hitInput: "nuslhitsCryoE"
+NuGraphCryoE.LoaderTool.spsInput: "nuslhitsCryoE"
+NuGraphCryoE.DecoderTools.SemanticDecoderTool.hitInput: "nuslhitsCryoE"
+NuGraphCryoE.avgs_u: [168.04594 , 178.3245  , 266.6149 , 3.857218 ]
+NuGraphCryoE.avgs_v: [1245.3547 , 176.54117 , 323.52786, 4.3267984]
+NuGraphCryoE.avgs_y: [1225.5012 , 183.58075 , 310.83493, 4.3409133]
+NuGraphCryoE.devs_u: [ 82.80644 ,  67.60649 , 274.32666, 1.2912455]
+NuGraphCryoE.devs_v: [ 293.06314,  66.8194  , 322.11386, 1.4249923]
+NuGraphCryoE.devs_y: [ 307.1943 ,  67.063324, 312.461  , 1.4532351]
+NuGraphCryoE.modelFileName: "model_mpvmpr_bnb_numu_cos.pt"
+#NuGraphCryoE.debug: true
+
+NuGraphCryoW: @local::NuGraphCryoE
+NuGraphCryoW.LoaderTool.hitInput: "nuslhitsCryoW"
+NuGraphCryoW.LoaderTool.spsInput: "nuslhitsCryoW"
+NuGraphCryoW.DecoderTools.SemanticDecoderTool.hitInput: "nuslhitsCryoW"
+
+END_PROLOG
diff --git a/icaruscode/TPC/NuGraph/scripts/CMakeLists.txt b/icaruscode/TPC/NuGraph/scripts/CMakeLists.txt
new file mode 100644
index 000000000..e0f7a8d2b
--- /dev/null
+++ b/icaruscode/TPC/NuGraph/scripts/CMakeLists.txt
@@ -0,0 +1,4 @@
+install_scripts(
+  "setup_tritonserver-nugraph-v0.sh"
+  "setup_tritonserver-nugraph-v0_grid.sh"
+  )
diff --git a/icaruscode/TPC/NuGraph/scripts/setup_tritonserver-nugraph-v0.sh b/icaruscode/TPC/NuGraph/scripts/setup_tritonserver-nugraph-v0.sh
new file mode 100755
index 000000000..31506929f
--- /dev/null
+++ b/icaruscode/TPC/NuGraph/scripts/setup_tritonserver-nugraph-v0.sh
@@ -0,0 +1,8 @@
+#!/bin/bash
+
+unset PYTHONHOME
+unset PYTHONPATH
+rm -f tritonserver_nugraph-v0.log
+export APPTAINER_BIND=/etc/hosts,/tmp,/cvmfs
+/cvmfs/oasis.opensciencegrid.org/mis/apptainer/1.3.2/bin/apptainer run --pid --ipc --home ~/:${HOME} --pwd ${PWD} /cvmfs/icarus.opensciencegrid.org/containers/tritonserver/nugraph-v0/ >& tritonserver_nugraph-v0.log &
+grep -q "Started" <(tail -f tritonserver_nugraph-v0.log)
diff --git a/icaruscode/TPC/NuGraph/scripts/setup_tritonserver-nugraph-v0_grid.sh b/icaruscode/TPC/NuGraph/scripts/setup_tritonserver-nugraph-v0_grid.sh
new file mode 100755
index 000000000..174a461d6
--- /dev/null
+++ b/icaruscode/TPC/NuGraph/scripts/setup_tritonserver-nugraph-v0_grid.sh
@@ -0,0 +1,31 @@
+#!/bin/bash
+#
+unset PYTHONHOME
+unset PYTHONPATH
+rm -f tritonserver_nugraph-v0.log
+export APPTAINER_BIND=/etc/hosts,/tmp,/cvmfs
+#
+export BASEPORT=8000
+export NPORTS=3
+#
+export HTTPPORT=$BASEPORT
+export GRPCPORT=$((BASEPORT+1))
+export METRPORT=$((BASEPORT+2))
+#
+if [ -z "$FCL" ]; then 
+    export FCL=$1
+fi
+#
+while 2>/dev/null >"/dev/tcp/0.0.0.0/$BASEPORT"; do
+    BASEPORT=$((BASEPORT+NPORTS))
+    export HTTPPORT=$BASEPORT
+    export GRPCPORT=$((BASEPORT+1))
+    export METRPORT=$((BASEPORT+2))
+done
+#
+echo "physics.producers.NuGraphCryoE.TritonConfig.serverURL: 'localhost:${GRPCPORT}'" >> "$FCL"
+echo "physics.producers.NuGraphCryoW.TritonConfig.serverURL: 'localhost:${GRPCPORT}'" >> "$FCL"
+#
+/cvmfs/oasis.opensciencegrid.org/mis/apptainer/1.3.2/bin/apptainer exec --pid --ipc --home ~/:${HOME} --pwd ${PWD} /cvmfs/icarus.opensciencegrid.org/containers/tritonserver/nugraph-v0/ tritonserver --model-repository /triton-server-config/models --http-port=${HTTPPORT} --grpc-port=${GRPCPORT} --metrics-port=${METRPORT} >& tritonserver_nugraph-v0.log &
+#
+grep -q "Started" <(tail -f tritonserver_nugraph-v0.log)
diff --git a/icaruscode/TPC/NuGraph/testinference_slice_icarus.fcl b/icaruscode/TPC/NuGraph/testinference_slice_icarus.fcl
new file mode 100644
index 000000000..84ceca649
--- /dev/null
+++ b/icaruscode/TPC/NuGraph/testinference_slice_icarus.fcl
@@ -0,0 +1,114 @@
+###
+### Test fcl file for NuGraph inference. For each cryostat it runs two producers (ICARUSNuSliceHitsProducer and NuGraphInference) plus one analyzer (NuGraphAnalyzer).
+### Author: G. Cerati (FNAL)
+###
+
+#include "nugraph_icarus.fcl"
+#include "services_common_icarus.fcl"
+#include "services_icarus_simulation.fcl"
+
+process_name: testinference
+
+services:
+{
+  TFileService:          { fileName: "reco_hist.root" }
+  TimeTracker:           {}
+  MemoryTracker:         {}
+  RandomNumberGenerator: {}
+  @table::icarus_common_services
+  @table::icarus_backtracking_services 
+  message:               @local::icarus_message_services_prod_debug
+  ParticleInventoryService: { 
+        ParticleInventory: {
+        EveIdCalculator: "EmEveIdCalculator"
+        G4ModuleLabel: "largeant"
+	OverrideRealData: true
+     }
+  }
+}
+
+source:
+{
+  module_type: RootInput
+  maxEvents:   -1
+}
+
+outputs:
+{
+ rootOutput:
+ {
+   module_type: RootOutput
+   dataTier: "reconstructed"
+   compressionLevel: 1
+   saveMemoryObjectThreshold: 0
+   fileName: "%ifb_%tc-%p.root"
+   fileProperties: {maxInputFiles: 1}
+   checkFileName: false
+   SelectEvents: []
+   outputCommands: [
+         "keep *_*_*_*"
+         #"drop *_*_*_*",
+         #"keep *_*_*_testinference"
+     ]
+ }
+}
+
+physics:
+{
+
+  producers:
+  {
+    nuslhitsCryoE: @local::nuslhitsCryoE
+    nuslhitsCryoW: @local::nuslhitsCryoW
+    #nuslhits: {
+    #  module_type:             "HitMerger"
+    #  HitProducerLabelVec:     ["nuslhitsCryoE","nuslhitsCryoW"]
+    #}
+    NuGraphCryoE: @local::NuGraphCryoE
+    NuGraphCryoW: @local::NuGraphCryoW
+  }
+  analyzers:
+  {
+    NuGraphAnaCryoE: { 
+      module_type: "NuGraphAnalyzer"
+      NuGraphLabel: "NuGraphCryoE"
+    }
+    NuGraphAnaCryoW: { 
+      module_type: "NuGraphAnalyzer"
+      NuGraphLabel: "NuGraphCryoW"
+    }
+  }
+
+  reco: [nuslhitsCryoE, nuslhitsCryoW, NuGraphCryoE, NuGraphCryoW]
+  #reco: [ nuslhitsCryoE, NuGraphCryoE]
+  trigger_paths:       [ reco ]
+
+  ana: [NuGraphAnaCryoE,NuGraphAnaCryoW]
+  #ana: [NuGraphAnaCryoE]
+  streamROOT:         [ rootOutput ]
+  end_paths:          [ ana, streamROOT ]
+}
+
+#physics.producers.nuslhitsCryoE.module_type: "ICARUSTrueNuSliceHitsProducer"
+#physics.producers.nuslhitsCryoW.module_type: "ICARUSTrueNuSliceHitsProducer"
+
+#physics.producers.NuGraphCryoE.debug: true
+#physics.producers.NuGraphCryoW.debug: true
+
+services.BackTrackerService.BackTracker.G4ModuleLabel: "largeant"
+services.BackTrackerService.BackTracker.SimChannelModuleLabel: "daq:simpleSC"
+
+
+### Parameters and instructions to test the NuSonic inference model (still under development)
+# in order to setup the triton server locally, please do:
+# > mrbslp
+# > setup_tritonservier-nugraph-v0.sh
+# then after the job you can kill the server process with:
+# > killall -9 /cvmfs/oasis.opensciencegrid.org/mis/apptainer/1.3.2/x86_64/libexec/apptainer/libexec/starter
+#NuGraphCryoE: @local::ApptainerNuGraphNuSonicTriton
+#NuGraphCryoE.LoaderTool.tool_type: "ICARUSNuGraphLoader"
+#NuGraphCryoE.LoaderTool.hitInput: "nuslhitsCryoE"
+#NuGraphCryoE.LoaderTool.spsInput: "nuslhitsCryoE"
+#NuGraphCryoE.DecoderTools.SemanticDecoderTool.hitInput: "nuslhitsCryoE"
+#NuGraphCryoE.TritonConfig.modelName: "nugraph2_icarus_mpvmprbnb"
+