diff --git a/fcl/reco/Definitions/stage1_icarus_defs.fcl b/fcl/reco/Definitions/stage1_icarus_defs.fcl
index dd82375cc..af83cde3d 100644
--- a/fcl/reco/Definitions/stage1_icarus_defs.fcl
+++ b/fcl/reco/Definitions/stage1_icarus_defs.fcl
@@ -16,8 +16,11 @@
 #include "flashmatch_simple_icarus.fcl"
 #include "crttrackproducer_icarus.fcl"
 #include "crtbacktracker_icarus.fcl"
+### FPoppi : The following CRT matching algorithms are deprecated.
 #include "crtt0matchingalg_icarus.fcl"
 #include "crtt0producer_icarus.fcl"
+### FPoppi : The following one is mantained.
+#include "crtt0tagging.fcl"
 ## The below can be found from the softlink to Supera in sbncode
 #include "supera_modules.fcl"
 #include "crtpmtmatching_parameters.fcl"
@@ -65,9 +68,11 @@ icarus_stage1_producers:
 
  ## crt producer
   crttrack:          @local::standard_crttrackproducer
-  CRTT0Matching:     @local::standard_crtt0producer
-  CRTT0MatchingW:    @local::standard_crtt0producerW
-  CRTT0MatchingE:    @local::standard_crtt0producerE
+  ### FPoppi: the following are deprecated (and the one above not mantained)
+  ## CRTT0Matching:     @local::standard_crtt0producer
+  ## CRTT0MatchingW:    @local::standard_crtt0producerW
+  ## CRTT0MatchingE:    @local::standard_crtt0producerE
+  CRTT0Tagging:         @local::icarus_crtt0tagging_data
 
   tpcpmtbarycentermatchCryoE: @local::data_tpcpmtbarycentermatchproducer_east
   tpcpmtbarycentermatchCryoW: @local::data_tpcpmtbarycentermatchproducer_west
@@ -223,6 +228,8 @@ icarus_crtt0match:                 [CRTT0Matching]
 
 icarus_crtt0match_eff:             [CRTT0MatchingE, CRTT0MatchingW]
 
+icarus_crtt0tagging:               [CRTT0Tagging]
+
 ### Below we include overrides for the modules above
 
 ## Overrides for filtering of cluster3D hits
diff --git a/fcl/reco/Stage1/Run1/stage1_multiTPC_icarus_gauss.fcl b/fcl/reco/Stage1/Run1/stage1_multiTPC_icarus_gauss.fcl
index d9ad45958..f8518077f 100644
--- a/fcl/reco/Stage1/Run1/stage1_multiTPC_icarus_gauss.fcl
+++ b/fcl/reco/Stage1/Run1/stage1_multiTPC_icarus_gauss.fcl
@@ -6,7 +6,7 @@ physics.reco: [ @sequence::icarus_filter_cluster3D,
 		@sequence::icarus_pandora_Gauss,
 		@sequence::icarus_reco_fm,
                 @sequence::icarus_crttrack,
-		@sequence::icarus_crtt0match,
+		@sequence::icarus_crtt0tagging,
                 caloskimCalorimetryCryoE, caloskimCalorimetryCryoW]
 
 physics.outana:            [ caloskimE, caloskimW, simpleLightAna, CRTDataAnalysis]
diff --git a/fcl/reco/Stage1/Run1/stage1_multiTPC_nofilter_icarus_gauss_crtpmt.fcl b/fcl/reco/Stage1/Run1/stage1_multiTPC_nofilter_icarus_gauss_crtpmt.fcl
index 39d2a83e2..2ac9bd5fb 100644
--- a/fcl/reco/Stage1/Run1/stage1_multiTPC_nofilter_icarus_gauss_crtpmt.fcl
+++ b/fcl/reco/Stage1/Run1/stage1_multiTPC_nofilter_icarus_gauss_crtpmt.fcl
@@ -6,7 +6,7 @@ physics.reco: [ @sequence::icarus_filter_cluster3D,
                 @sequence::icarus_pandora_Gauss,
 		@sequence::icarus_reco_fm,
                 @sequence::icarus_crttrack,
-		@sequence::icarus_crtt0match,
+		@sequence::icarus_crtt0tagging,
                 caloskimCalorimetryCryoE, caloskimCalorimetryCryoW]
 
 physics.outana:            [ caloskimE, caloskimW, simpleLightAna, CRTDataAnalysis, CRTAnalysis]
diff --git a/fcl/reco/Stage1/Run2/stage1_run2_icarus.fcl b/fcl/reco/Stage1/Run2/stage1_run2_icarus.fcl
index 40ea40c2a..fdd161499 100644
--- a/fcl/reco/Stage1/Run2/stage1_run2_icarus.fcl
+++ b/fcl/reco/Stage1/Run2/stage1_run2_icarus.fcl
@@ -7,7 +7,7 @@ physics.reco: [ @sequence::icarus_filter_cluster3D,
                 @sequence::icarus_reco_fm,
                 @sequence::icarus_tpcpmtbarycentermatch,
                 @sequence::icarus_crttrack,
-                @sequence::icarus_crtt0match,
+                @sequence::icarus_crtt0tagging,
                 caloskimCalorimetryCryoE, caloskimCalorimetryCryoW]
 
 physics.outana:            [ @sequence::icarus_analysis_modules ]
diff --git a/icaruscode/CRT/CMakeLists.txt b/icaruscode/CRT/CMakeLists.txt
index 8bfd2d749..5c6e34abf 100644
--- a/icaruscode/CRT/CMakeLists.txt
+++ b/icaruscode/CRT/CMakeLists.txt
@@ -1,7 +1,7 @@
 
 add_subdirectory(CRTUtils)
 add_subdirectory(CRTDecoder)
-
+add_subdirectory(CRTLegacyCode)
 
 art_make(
     NO_PLUGINS
@@ -11,7 +11,6 @@ art_make(
         CRTDetSim_module.cc
         CRTSimHitProducer_module.cc
         CRTTrueHitProducer_module.cc
-        CRTTzeroProducer_module.cc
         CRTTrackProducer_module.cc
         CRTSimAnalysis_module.cc
         CRTDataAnalysis_module.cc
@@ -20,7 +19,6 @@ art_make(
         CRTAutoVeto_module.cc
         FlashResAna_module.cc
         PhotBackground_module.cc
-        CRTT0Matching_module.cc
         CRTTPCMatchingAna_module.cc	
         CRTPMTMatchingAna_module.cc
         CRTTPCTruthEff_module.cc
@@ -62,6 +60,25 @@ art_make_library(
                 cetlib::cetlib
         )
 
+art_make_library(
+        LIBRARY_NAME
+                icaruscode_CRTMatchingUtils
+        SOURCE
+                CRTUtils/CRTMatchingUtils.cxx
+        LIBRARIES
+                fhiclcpp::fhiclcpp
+                canvas::canvas
+                lardataobj::RecoBase
+                lardata::Utilities
+	        lardata::ArtDataHelper
+                larcorealg::Geometry
+                lardataalg::DetectorInfo
+                sbnobj::Common_CRT
+                cetlib::cetlib
+                messagefacility::MF_MessageLogger
+                ROOT::Matrix
+        )
+
 cet_build_plugin(CRTGeometryHelper art::service
      LIBRARIES
         larcorealg::Geometry
@@ -129,23 +146,6 @@ cet_build_plugin( CRTTrueHitProducer art::module
         cetlib::cetlib
         )
 
-cet_build_plugin(CRTTzeroProducer art::module
-     LIBRARIES
-        larcorealg::Geometry
-        icaruscode_CRT
-        sbnobj::ICARUS_CRT
-        sbnobj::Common_CRT
-        art_root_io::TFileService_service
-        lardataalg::DetectorInfo
-        nurandom::RandomUtils_NuRandomService_service
-        art::Framework_Services_Registry
-        art::Framework_Services_Optional_RandomNumberGenerator_service
-        messagefacility::MF_MessageLogger
-        messagefacility::headers
-        CLHEP::CLHEP
-        lardata::Utilities
-        )
-
 cet_build_plugin(CRTTrackProducer art::module
      LIBRARIES
         larcorealg::Geometry
@@ -382,78 +382,6 @@ cet_build_plugin( PhotBackground art::module
         ROOT::Tree
         )
 
-cet_build_plugin(CRTT0Matching art::module
-     LIBRARIES
-              icaruscode_CRTData
-              icaruscode_CRT
-              sbnobj::Common_CRT
-              icaruscode_CRTUtils
-              larcore::Geometry_Geometry_service
-              larsim::Simulation lardataobj::Simulation
-              larsim::MCCheater_BackTrackerService_service
-              larsim::MCCheater_ParticleInventoryService_service
-              lardata::Utilities
-              larevt::Filters
-              lardataobj::RawData
-              lardataobj::RecoBase
-              lardataobj::AnalysisBase
-              lardata::RecoObjects
-              larpandora::LArPandoraInterface
-              larcorealg::Geometry
-              nusimdata::SimulationBase
-              art::Framework_Core
-              art::Framework_Principal
-              art::Framework_Services_Registry
-              art_root_io::tfile_support
-              art_root_io::TFileService_service
-              art::Persistency_Common canvas::canvas
-              art::Persistency_Provenance
-              art::Utilities
-              messagefacility::MF_MessageLogger
-              ROOT::Core
-              ROOT::Geom
-              ROOT::XMLIO
-              ROOT::Gdml
-              ROOT::Tree
-              ROOT::Spectrum
-              ROOT::RooFit
-              ROOT::RooFitCore
-)
-
-
-cet_build_plugin(CRTT0MatchingAna art::module
-     LIBRARIES
-              sbnobj::ICARUS_CRT
-              icaruscode_CRT
-              sbnobj::Common_CRT
-              icaruscode_CRTUtils
-              larcorealg::Geometry
-              larcore::Geometry_Geometry_service
-              larsim::Simulation lardataobj::Simulation
-              larsim::MCCheater_BackTrackerService_service
-              larsim::MCCheater_ParticleInventoryService_service
-              lardata::Utilities
-              larevt::Filters
-              lardataobj::RawData
-              lardataobj::RecoBase
-              lardataobj::AnalysisBase
-              lardata::RecoObjects
-              larpandora::LArPandoraInterface
-              larcorealg::Geometry
-              nusimdata::SimulationBase
-              art::Persistency_Common canvas::canvas
-              art::Persistency_Provenance
-              art::Utilities
-              ROOT::Core
-              ROOT::Tree
-              ROOT::Geom
-              ROOT::XMLIO
-              ROOT::Gdml
-              ROOT::Spectrum
-              ROOT::RooFit
-              ROOT::RooFitCore
-)
-
 cet_build_plugin(CRTTPCMatchingAna art::module
      LIBRARIES
               sbnobj::ICARUS_CRT
@@ -517,6 +445,26 @@ simple_plugin(CRTTPCTruthEff module
 
 )
 
+cet_build_plugin(CRTT0Tagging art::module
+     LIBRARIES
+        icaruscode::IcarusObj
+        icaruscode_CRTMatchingUtils
+        icaruscode_CRTUtils
+        sbnobj::Common_CRT
+        larcorealg::Geometry
+        lardataalg::DetectorInfo
+        larcore::Geometry_Geometry_service
+        larsim::Simulation lardataobj::Simulation
+        larsim::MCCheater_BackTrackerService_service
+        larsim::MCCheater_ParticleInventoryService_service
+        lardata::headers
+        lardataobj::RecoBase
+        lardataobj::AnalysisBase
+        nusimdata::SimulationBase
+        art_root_io::TFileService_service
+        art_root_io::tfile_support
+        ROOT::Tree
+)
 
 install_headers()
 install_fhicl()
diff --git a/icaruscode/CRT/CRTLegacyCode/CMakeLists.txt b/icaruscode/CRT/CRTLegacyCode/CMakeLists.txt
new file mode 100644
index 000000000..3d6fda488
--- /dev/null
+++ b/icaruscode/CRT/CRTLegacyCode/CMakeLists.txt
@@ -0,0 +1,93 @@
+cet_build_plugin(CRTTzeroProducer art::module
+     LIBRARIES
+        larcorealg::Geometry
+        icaruscode_CRT
+        sbnobj::ICARUS_CRT
+        sbnobj::Common_CRT
+        art_root_io::TFileService_service
+        lardataalg::DetectorInfo
+        nurandom::RandomUtils_NuRandomService_service
+        art::Framework_Services_Registry
+        art::Framework_Services_Optional_RandomNumberGenerator_service
+        messagefacility::MF_MessageLogger
+        messagefacility::headers
+        CLHEP::CLHEP
+        lardata::Utilities
+        )
+
+cet_build_plugin(CRTT0Matching art::module
+     LIBRARIES
+              icaruscode_CRTData
+              icaruscode_CRT
+              sbnobj::Common_CRT
+              icaruscode_CRTUtils
+              larcore::Geometry_Geometry_service
+              larsim::Simulation lardataobj::Simulation
+              larsim::MCCheater_BackTrackerService_service
+              larsim::MCCheater_ParticleInventoryService_service
+              lardata::Utilities
+              larevt::Filters
+              lardataobj::RawData
+              lardataobj::RecoBase
+              lardataobj::AnalysisBase
+              lardata::RecoObjects
+              larpandora::LArPandoraInterface
+              larcorealg::Geometry
+              nusimdata::SimulationBase
+              art::Framework_Core
+              art::Framework_Principal
+              art::Framework_Services_Registry
+              art_root_io::tfile_support
+              art_root_io::TFileService_service
+              art::Persistency_Common canvas::canvas
+              art::Persistency_Provenance
+              art::Utilities
+              messagefacility::MF_MessageLogger
+              ROOT::Core
+              ROOT::Geom
+              ROOT::XMLIO
+              ROOT::Gdml
+              ROOT::Tree
+              ROOT::Spectrum
+              ROOT::RooFit
+              ROOT::RooFitCore
+)
+
+
+cet_build_plugin(CRTT0MatchingAna art::module
+     LIBRARIES
+              sbnobj::ICARUS_CRT
+              icaruscode_CRT
+              sbnobj::Common_CRT
+              icaruscode_CRTUtils
+              larcorealg::Geometry
+              larcore::Geometry_Geometry_service
+              larsim::Simulation lardataobj::Simulation
+              larsim::MCCheater_BackTrackerService_service
+              larsim::MCCheater_ParticleInventoryService_service
+              lardata::Utilities
+              larevt::Filters
+              lardataobj::RawData
+              lardataobj::RecoBase
+              lardataobj::AnalysisBase
+              lardata::RecoObjects
+              larpandora::LArPandoraInterface
+              larcorealg::Geometry
+              nusimdata::SimulationBase
+              art::Persistency_Common canvas::canvas
+              art::Persistency_Provenance
+              art::Utilities
+              ROOT::Core
+              ROOT::Tree
+              ROOT::Geom
+              ROOT::XMLIO
+              ROOT::Gdml
+              ROOT::Spectrum
+              ROOT::RooFit
+              ROOT::RooFitCore
+)
+
+
+install_headers()
+install_fhicl()
+install_source()
diff --git a/icaruscode/CRT/CRTT0MatchingAna_module.cc b/icaruscode/CRT/CRTLegacyCode/CRTT0MatchingAna_module.cc
similarity index 100%
rename from icaruscode/CRT/CRTT0MatchingAna_module.cc
rename to icaruscode/CRT/CRTLegacyCode/CRTT0MatchingAna_module.cc
diff --git a/icaruscode/CRT/CRTT0Matching_module.cc b/icaruscode/CRT/CRTLegacyCode/CRTT0Matching_module.cc
similarity index 100%
rename from icaruscode/CRT/CRTT0Matching_module.cc
rename to icaruscode/CRT/CRTLegacyCode/CRTT0Matching_module.cc
diff --git a/icaruscode/CRT/CRTTzeroProducer_module.cc b/icaruscode/CRT/CRTLegacyCode/CRTTzeroProducer_module.cc
similarity index 100%
rename from icaruscode/CRT/CRTTzeroProducer_module.cc
rename to icaruscode/CRT/CRTLegacyCode/CRTTzeroProducer_module.cc
diff --git a/icaruscode/CRT/crtt0matching_icarus.fcl b/icaruscode/CRT/CRTLegacyCode/crtt0matching_icarus.fcl
similarity index 100%
rename from icaruscode/CRT/crtt0matching_icarus.fcl
rename to icaruscode/CRT/CRTLegacyCode/crtt0matching_icarus.fcl
diff --git a/icaruscode/CRT/crtt0matchingana_icarus.fcl b/icaruscode/CRT/CRTLegacyCode/crtt0matchingana_icarus.fcl
similarity index 100%
rename from icaruscode/CRT/crtt0matchingana_icarus.fcl
rename to icaruscode/CRT/CRTLegacyCode/crtt0matchingana_icarus.fcl
diff --git a/icaruscode/CRT/crtt0producer_icarus.fcl b/icaruscode/CRT/CRTLegacyCode/crtt0producer_icarus.fcl
similarity index 100%
rename from icaruscode/CRT/crtt0producer_icarus.fcl
rename to icaruscode/CRT/CRTLegacyCode/crtt0producer_icarus.fcl
diff --git a/icaruscode/CRT/crttzeroproducer_icarus.fcl b/icaruscode/CRT/CRTLegacyCode/crttzeroproducer_icarus.fcl
similarity index 100%
rename from icaruscode/CRT/crttzeroproducer_icarus.fcl
rename to icaruscode/CRT/CRTLegacyCode/crttzeroproducer_icarus.fcl
diff --git a/icaruscode/CRT/CRTT0Tagging_module.cc b/icaruscode/CRT/CRTT0Tagging_module.cc
new file mode 100644
index 000000000..166dc4407
--- /dev/null
+++ b/icaruscode/CRT/CRTT0Tagging_module.cc
@@ -0,0 +1,589 @@
+/**
+ * @file   icaruscode/CRT/CRTT0Tagging_module.cc
+ * @author Francesco Poppi (poppi@bo.infn.it)
+ * @date   October 2024
+ */
+
+#include "sbnobj/Common/CRT/CRTHit.hh"
+#include "icaruscode/CRT/CRTUtils/CRTMatchingUtils.h"
+#include "sbnobj/Common/CRT/CRTHitT0TaggingInfo.hh"
+#include "sbnobj/Common/CRT/CRTHitT0TaggingTruthInfo.hh"
+#include "icaruscode/CRT/CRTUtils/RecoUtils.h"
+#include "icaruscode/CRT/CRTUtils/CRTCommonUtils.h"
+// Framework includes
+#include "art/Framework/Core/EDProducer.h"
+#include "art/Framework/Core/ModuleMacros.h"
+#include "art/Framework/Principal/Handle.h"
+#include "art/Framework/Principal/Event.h" 
+#include "art/Framework/Principal/Run.h"
+#include "canvas/Utilities/Exception.h"
+#include "canvas/Persistency/Common/Ptr.h" 
+#include "canvas/Persistency/Common/PtrVector.h" 
+#include "art/Persistency/Common/PtrMaker.h"
+#include "canvas/Persistency/Common/FindManyP.h"
+#include "canvas/Persistency/Common/FindOne.h"
+#include "art_root_io/TFileService.h"
+#include "canvas/Utilities/InputTag.h"
+#include "fhiclcpp/ParameterSet.h"
+#include "messagefacility/MessageLogger/MessageLogger.h"
+#include <memory>
+#include <vector>
+#include <map>
+#include <algorithm>
+#include <cmath> // std::abs(), std::hypot()
+// LArSoft
+#include "larcorealg/Geometry/GeometryCore.h"
+#include "larcorealg/CoreUtils/zip.h"
+#include "lardataobj/AnalysisBase/T0.h"
+#include "lardataobj/RecoBase/Hit.h"
+#include "lardataobj/RecoBase/Track.h"
+#include "lardataobj/RecoBase/PFParticleMetadata.h"
+#include "lardataobj/RecoBase/PFParticle.h"
+#include "larcore/CoreUtils/ServiceUtil.h" // lar::providerFrom()
+#include "larcore/Geometry/Geometry.h"
+#include "lardata/DetectorInfoServices/LArPropertiesService.h"
+#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
+#include "larsim/MCCheater/BackTrackerService.h"
+#include "larsim/MCCheater/ParticleInventoryService.h"
+// ROOT
+#include "TTree.h"
+
+namespace icarus::crt {
+
+class CRTT0Tagging;
+
+}  // namespace icarus::crt
+
+using namespace icarus::crt;
+
+class icarus::crt::CRTT0Tagging : public art::EDProducer {
+public:
+  
+  using CRTHit = sbn::crt::CRTHit;
+  //using CRTPMTMatching = sbn::crt::CRTPMTMatching;
+
+  explicit CRTT0Tagging(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.
+  CRTT0Tagging(CRTT0Tagging const&) = delete;
+  CRTT0Tagging(CRTT0Tagging&&) = delete;
+  CRTT0Tagging& operator=(CRTT0Tagging const&) = delete;
+  CRTT0Tagging& operator=(CRTT0Tagging&&) = delete;
+
+  // Required functions.
+  bool hasModID(std::uint32_t modID, sbn::crt::CRTHit const& crthit);
+  //void beginRun(art::Run& r) override;
+  void beginJob() override;
+  void produce(art::Event& e) override;
+
+private:
+
+  // Declare member data here.
+
+  // MC Truth
+  art::InputTag fSimulationProducerLabel; 
+  art::InputTag fAuxDetSimProducerLabel;
+  //art::InputTag fCRTSimHitProducerLabel; 
+  //art::InputTag fCRTTrueHitProducerLabel;
+  //art::InputTag fCRTDetSimProducerLabel;
+  //art::InputTag fCRTSimTrackProducerLabel;
+  art::InputTag fSimChannelProducerLabel;
+
+  art::InputTag fCrtHitModuleLabel;
+  
+  std::vector<art::InputTag> fTPCTrackLabel; ///< labels for source of tracks
+  std::vector<art::InputTag> fPFParticleLabel; ///< labels for source of PFParticle
+  std::vector<art::InputTag> fHitLabel; ///< labels for source of hits
+  std::vector<art::InputTag> fTRKHMLabel; ///< labels for hit metadata
+  
+  art::ServiceHandle<art::TFileService> tfs;
+
+  CRTCommonUtils fCrtUtils; 
+  CRTMatchingAlg fMatchingAlg;
+  geo::GeometryCore const* fGeometryService;  ///< pointer to Geometry provider
+  
+  double fMinimalTrackLength;
+  int fMinimumGoodHits;
+  double fMaximalCRTDistance;
+  double fGoodCandidateDistance;
+  double fMaximumDeltaT;
+  bool fData;
+  bool fSkipTruth;
+
+  icarus::crt::dataTools::TopCRTCentersMap fTopCRTCenterMap;
+  icarus::crt::dataTools::TopCRTTransformations fTopCRTTransformations;
+
+  TTree* fTree;
+  int fEvent;        ///< number of the event being processed
+  int fRun;          ///< number of the run being processed
+  int fSubRun;       ///< number of the sub-run being processed
+  int fCrtRegion;
+  int fCrtSys;
+  int fCryo;
+  double fTrackLength;
+  double fFirstEigenValue;
+  double fSecondEigenValue;
+  double fThirdEigenValue;
+  double fTrackCrtDistance;
+  double fTrackCrtDeltaX;
+  double fTrackCrtDeltaY;
+  double fTrackCrtDeltaZ;
+  double fCrtX;
+  double fCrtY;
+  double fCrtZ;
+  double fCrossPointX;
+  double fCrossPointY;
+  double fCrossPointZ;      
+  bool fTrueMatch;
+
+};
+
+icarus::crt::CRTT0Tagging::CRTT0Tagging(fhicl::ParameterSet const& p)
+    : EDProducer{p},
+      fSimulationProducerLabel(p.get<art::InputTag>("SimulationLabel","largeant")),
+      fAuxDetSimProducerLabel(p.get<art::InputTag>("AuxDetSimProducerLabel","genericcrt")),
+      fSimChannelProducerLabel(p.get<art::InputTag>("SimChannelProducer", {"daq:simpleSC"})),
+      fCrtHitModuleLabel(p.get<art::InputTag>("CrtHitModuleLabel", "crthit")),
+      fTPCTrackLabel(p.get< std::vector<art::InputTag> >("TPCTrackLabel",             {""})),
+      fPFParticleLabel(p.get< std::vector<art::InputTag> >("PFParticleLabel",             {""})),  
+      fHitLabel(p.get< std::vector<art::InputTag> >("HitLabel",             {""})),
+      fTRKHMLabel(p.get< std::vector<art::InputTag> > ("TRKHMLabel", {""})),
+      fMatchingAlg(p.get<fhicl::ParameterSet> ("MatchingAlg")),
+      fGeometryService(lar::providerFrom<geo::Geometry>()),
+      fMinimalTrackLength(p.get<double>("MinimalTrackLength", 40.0)),
+      fMinimumGoodHits(p.get<double>("MinimumGoodHits", 5)),
+      fMaximalCRTDistance(p.get<double>("MaximalCRTDistance", 300.)),
+      fGoodCandidateDistance(p.get<double>("GoodCandidateDistance", 100.)),
+      fMaximumDeltaT(p.get<double>("MaximumDeltaT", 10000.)),
+      fData(p.get<bool>("isData", true)),
+      fSkipTruth(p.get<bool>("skipTruth", false))
+{
+  
+  produces< std::vector<anab::T0>                   >();
+  produces< art::Assns<recob::Track , anab::T0>     >();
+  produces< art::Assns<sbn::crt::CRTHit, anab::T0>  >();  
+  produces< std::vector<sbn::crt::CRTHitT0TaggingInfo> >();
+  produces< art::Assns<recob::Track , sbn::crt::CRTHitT0TaggingInfo> >();
+  produces< art::Assns<recob::PFParticle , anab::T0>  >();
+  produces< art::Assns<recob::PFParticle , sbn::crt::CRTHitT0TaggingInfo> >();
+  produces< std::vector<sbn::crt::CRTHitT0TaggingTruthInfo> >();
+  produces< art::Assns<sbn::crt::CRTHitT0TaggingInfo, sbn::crt::CRTHitT0TaggingTruthInfo> >();
+
+  //produces< art::Assns<icarus::CRTTPCMatchingInfo, anab::T0>  >();  
+  //produces< art::Assns<recob::Track, icarus::CRTTPCMatchingInfo>  >();
+  //produces< art::Assns<sbn::crt::CRTHit, icarus::CRTTPCMatchingInfo>  >();
+
+  if (fTPCTrackLabel.size() != fPFParticleLabel.size()) {
+    throw art::Exception{ art::errors::Configuration }
+      << fTPCTrackLabel.size() << " TPC track data products configured (`TPCTrackLabel`), should have been "
+      << fPFParticleLabel.size();
+  }
+  if (fHitLabel.size() != fPFParticleLabel.size()) {
+    throw art::Exception{ art::errors::Configuration }
+      << fHitLabel.size() << " TPC hit data products configured (`HitLabel`), should have been "
+      << fPFParticleLabel.size();
+  }
+
+  if (fTRKHMLabel.size() > fPFParticleLabel.size()) {
+    throw art::Exception{ art::errors::Configuration }
+      << fTRKHMLabel.size() << " track-hit metadata data products configured (`TRKHMproducer`), should have been "
+      << fPFParticleLabel.size();
+  }
+  fTRKHMLabel.resize(fPFParticleLabel.size()); // extend with empty labels
+  // replace empty defaults with the actual input tag value, assumed the same as tracks
+  for (std::size_t i = 0; i < fTRKHMLabel.size(); ++i)
+    if (fTRKHMLabel[i].empty()) fTRKHMLabel[i] = fTPCTrackLabel[i];
+
+}
+
+bool icarus::crt::CRTT0Tagging::hasModID(std::uint32_t modID, sbn::crt::CRTHit const& crthit){
+  for(auto const& mactopes: crthit.pesmap){
+    for(auto const& chanpe: mactopes.second){
+      int thisModID=(int)fCrtUtils.MacToAuxDetID(mactopes.first, chanpe.first);
+      if(thisModID==(int)modID) return true;
+    }
+  }
+  return false;
+}
+
+//void icarus::crt::CRTT0Tagging::beginRun(art::Run& r)
+void icarus::crt::CRTT0Tagging::beginJob()
+{
+  fTopCRTCenterMap=icarus::crt::dataTools::LoadTopCRTCenters();
+  fTopCRTTransformations=icarus::crt::dataTools::LoadTopCRTTransformations();
+
+  fTree = tfs->make<TTree>("matchTree","CRTHit - TPC track matching analysis");
+
+  fTree->Branch("Event",              &fEvent,              "Event/I");
+  fTree->Branch("SubRun",             &fSubRun,             "SubRun/I");
+  fTree->Branch("Run",                &fRun,                "Run/I");
+  fTree->Branch("Cryo",               &fCryo,               "Cryo/I");
+  fTree->Branch("CrtSys",             &fCrtSys,             "CrtSys/I");
+  fTree->Branch("CrtRegion",          &fCrtRegion,          "CrtRegion/I");
+  fTree->Branch("TrackLength",        &fTrackLength        );
+  fTree->Branch("FirstEigenValue",    &fFirstEigenValue    );
+  fTree->Branch("SecondEigenValue",   &fSecondEigenValue   );
+  fTree->Branch("ThirdEigenValue",    &fThirdEigenValue    );
+  fTree->Branch("TrackCrtDistance",   &fTrackCrtDistance   );
+  fTree->Branch("TrackCrtDeltaX",     &fTrackCrtDeltaX     );
+  fTree->Branch("TrackCrtDeltaY",     &fTrackCrtDeltaY     );
+  fTree->Branch("TrackCrtDeltaZ",     &fTrackCrtDeltaZ     );
+  fTree->Branch("CrtX",               &fCrtX               );
+  fTree->Branch("CrtY",               &fCrtY               );
+  fTree->Branch("CrtZ",               &fCrtZ               );
+  fTree->Branch("CrossPointX",        &fCrossPointX        );
+  fTree->Branch("CrossPointY",        &fCrossPointY        );
+  fTree->Branch("CrossPointZ",        &fCrossPointZ        );
+  fTree->Branch("TrueMatch",          &fTrueMatch          );
+}
+
+void icarus::crt::CRTT0Tagging::produce(art::Event& e)
+{
+  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(e);
+  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(e, clockData);
+
+  mf::LogDebug("CRTT0Tagging: ") << "beginning production" << '\n';
+
+  auto t0col = std::make_unique< std::vector<anab::T0> > ();
+  auto trackAssn = std::make_unique< art::Assns<recob::Track, anab::T0> >();
+  auto t0CrtHitAssn = std::make_unique< art::Assns<sbn::crt::CRTHit, anab::T0> >();
+  auto matchInfoCol = std::make_unique< std::vector<sbn::crt::CRTHitT0TaggingInfo> >();
+  auto matchInfoTruthCol = std::make_unique< std::vector<sbn::crt::CRTHitT0TaggingTruthInfo> >();
+
+  auto pfpAssn = std::make_unique< art::Assns<recob::PFParticle, anab::T0> >();
+  auto trackMatchInfoAssn = std::make_unique< art::Assns<recob::Track , sbn::crt::CRTHitT0TaggingInfo> >();
+  auto pfpMatchInfoAssn = std::make_unique< art::Assns<recob::PFParticle , sbn::crt::CRTHitT0TaggingInfo> >();
+
+  auto truthAssn = std::make_unique< art::Assns<sbn::crt::CRTHitT0TaggingInfo, sbn::crt::CRTHitT0TaggingTruthInfo> >();
+
+  art::PtrMaker<anab::T0> makeT0ptr{ e }; // create art pointers to the new T0 
+  art::PtrMaker<sbn::crt::CRTHitT0TaggingInfo> makeMatchInfoPtr{ e }; // create art pointers to the CRTHitT0TaggingInfo
+  art::PtrMaker<sbn::crt::CRTHitT0TaggingTruthInfo> makeMatchTruthInfoPtr{ e }; // create art pointers to the CRTHitT0TaggingTruthInfo
+
+  std::map< int, const simb::MCParticle*> particleMap;
+  std::map<std::pair<int,double>,std::vector<int>> crtParticleMap;
+  std::map<int,bool> isNuMap;
+  std::vector<art::Ptr<sim::SimChannel>> simchannels;
+
+  // CRTHits
+  std::vector<art::Ptr<CRTHit>> CRTHitList;
+  art::ValidHandle<std::vector<CRTHit>> crthits = e.getValidHandle<std::vector<CRTHit>>(fCrtHitModuleLabel);
+  art::fill_ptr_vector(CRTHitList, crthits);
+
+  // If it is not data is MC.
+  // Retrieving MC truth information.
+  // Three maps (particleMap, crtParticleMap and isNuMap) are filled.
+  // -> particleMap maps truth level informations for all the MC particles. The key is the Geant4 TrackID
+  // -> crtParticleMap maps the CRTHits at truth level. The key is a pair of reconstructed CRT Hit module and time.
+  // The object of the map is a vector of Geant4 Track IDs of the particles that generated that hit.
+  // -> isNuMap maps the reconstructed tracks. The key is the Track ID and the object is a boolean True/False if
+  // the Track was a neutrino related interaction or not.
+  if(!fData){
+    if(fSkipTruth){
+      mf::LogInfo("CRTT0Tagging") <<"This is MC, but MC truth is not considered!";
+    } else{
+      mf::LogInfo("CRTT0Tagging") <<"This is MC, MC truth is considered!";
+    }   
+  }
+  if(!fData && !fSkipTruth){
+    art::ServiceHandle<cheat::ParticleInventoryService> partInventory;
+
+    art::Handle< std::vector<sim::SimChannel>> simChannelHandle;
+    if (!e.getByLabel(fSimChannelProducerLabel, simChannelHandle)){
+	    throw cet::exception("CRTT0Tagging") 
+	        << " No sim::SimChannel objects in this event - "
+	        << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
+    }
+    art::fill_ptr_vector(simchannels, simChannelHandle);
+
+    // Define "handle" to Generator level MCTruth objects
+    art::Handle< vector<simb::MCTruth>> genHandle;
+    // Define a "handle" to point to a vector of MCParticle objects.
+    art::Handle< vector<simb::MCParticle> > particleHandle;
+
+    if (!e.getByLabel("generator", genHandle)) {
+        std::cout << "could not get handle to gen objects!!!" << std::endl;
+    }
+
+    if (!e.getByLabel(fSimulationProducerLabel, particleHandle)) {
+	    // If we have no MCParticles at all in an event, but we are requiring
+      // to have MC truth information, throw exception.
+	    throw cet::exception("CRTT0Tagging") 
+	      << " No simb::MCParticle objects in this event - "
+	      << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
+    }
+
+    // Handle to AuxDetSimChannel (CRT module) objects generated by LArG4
+    art::Handle<vector<sim::AuxDetSimChannel> > auxDetSimChannelHandle;
+    if (!e.getByLabel(fAuxDetSimProducerLabel, auxDetSimChannelHandle)) {
+      throw cet::exception("CRTT0Tagging")
+        << " No sim::AuxDetSimChannel objects in this event - "
+        << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
+    }
+    //if((*genHandle).size()>1) 
+    //      throw cet::exception("CRTT0Tagging") << "gen stage MCParticle vector has more than 1 entry!" << '\n';
+    for ( auto const& particle : (*particleHandle) ){
+      // Add the address of the MCParticle to the map, with the
+      // track ID as the key.
+      particleMap[particle.TrackId()] = &particle;
+      art::Ptr<simb::MCTruth> mcTruth=partInventory->ParticleToMCTruth_P(&particle);
+      
+      bool isNu=false;
+      
+      if (mcTruth->Origin() == simb::kBeamNeutrino) isNu=true;
+
+      isNuMap[particle.TrackId()] = isNu;
+
+      for ( auto const& channel : (*auxDetSimChannelHandle) ){
+        auto const& auxDetIDEs = channel.AuxDetIDEs();
+        for ( auto const& ide : auxDetIDEs ){
+          if ( ide.trackID != particle.TrackId() ) continue;
+          if ( ide.energyDeposited * 1.0e6 < 50 ) continue; // skip energy deposits of less then 50 keV
+          size_t adid = channel.AuxDetID();
+          uint32_t region=fCrtUtils.AuxDetRegionNameToNum(fCrtUtils.GetAuxDetRegion(adid));
+          uint32_t modID=channel.AuxDetID();
+          float aveT = (ide.entryT + ide.exitT) / 2.0;
+          for(auto const& crthit : CRTHitList){
+            if(crthit->plane!=(int)region) continue;
+            if(abs(aveT-crthit->ts1_ns)>200) continue;
+            bool modFound=hasModID(modID, *crthit);
+            if(!modFound)continue;
+            std::pair<int,double> thisMatch = std::make_pair((int)crthit->feb_id[0],crthit->ts1_ns);
+            crtParticleMap[thisMatch].push_back(particle.TrackId());
+          } // CRT Hits loop
+        } // Energy deposits loop
+      } // CRT sim channels loop
+    } // MC particles loop
+  } // End MC Only
+  for(const auto& [ PFPLabel, TPCTrackLabel, HitLabel, TRKHMLabel ]: util::zip(fPFParticleLabel, fTPCTrackLabel, fHitLabel, fTRKHMLabel)){
+    std::vector<art::Ptr<recob::PFParticle>> PFParticleList;
+    art::ValidHandle<std::vector<recob::PFParticle>> pfparticles = e.getValidHandle<std::vector<recob::PFParticle>>(PFPLabel);
+    art::fill_ptr_vector(PFParticleList, pfparticles);
+
+    // Pandora MetaData
+    art::FindOne<anab::T0> fmt0pandora(pfparticles, e, PFPLabel);
+    art::FindOne<larpandoraobj::PFParticleMetadata> PFPMetaDataAssoc(pfparticles, e, PFPLabel);
+
+    // Tracks
+    art::ValidHandle<std::vector<recob::Track>> tracks = e.getValidHandle<std::vector<recob::Track>>(TPCTrackLabel);
+
+    // Track - associated data
+    art::FindManyP<recob::Track> fmTracks(PFParticleList, e, TPCTrackLabel);
+
+    // Collect all hits
+    art::ValidHandle<std::vector<recob::Hit>> allhit_handle = e.getValidHandle<std::vector<recob::Hit>>(HitLabel);
+    std::vector<art::Ptr<recob::Hit>> allHits;
+    art::fill_ptr_vector(allHits, allhit_handle);
+
+    std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> id_to_ide_map;
+    std::map<int, std::vector<art::Ptr<recob::Hit>>> id_to_truehit_map;
+    if(!simchannels.empty() && !fData){
+      art::ServiceHandle<cheat::BackTrackerService> btServ;
+      // ID (TrackID) refers to the reconstructed TrackID, it is not the Geant4 ID.
+      id_to_ide_map = RecoUtils::PrepSimChannels(simchannels, *fGeometryService);
+      id_to_truehit_map = RecoUtils::buildTrackIDtoHitsMap(allHits, clockData, *btServ.get());
+    }
+    // Start looping on the particles
+    for (art::Ptr<recob::PFParticle> const& p_pfp: PFParticleList) {
+      const std::vector<art::Ptr<recob::Track>> thisTrack = fmTracks.at(p_pfp.key());
+      if (thisTrack.size() != 1) continue;    
+      art::Ptr<recob::Track> trkPtr = thisTrack.at(0);
+      const recob::Track &track = *trkPtr;
+      if(track.Length()<fMinimalTrackLength) continue;
+      art::InputTag thm_label = TRKHMLabel;
+      art::FindManyP<recob::Hit, recob::TrackHitMeta> fmtrkHits(tracks, e, thm_label);
+      std::vector<art::Ptr<recob::Hit>> emptyHitVector;
+      const std::vector<art::Ptr<recob::Hit>> &trkHits  = fmtrkHits.isValid() ? fmtrkHits.at(trkPtr.key()) : emptyHitVector;
+      std::vector<const recob::TrackHitMeta*> emptyTHMVector;
+      const std::vector<const recob::TrackHitMeta*> &trkHitMetas = fmtrkHits.isValid() ? fmtrkHits.data(trkPtr.key()) : emptyTHMVector;
+      int trueTrackId=-9;
+      if(!fData) trueTrackId= abs(RecoUtils::TrueParticleIDFromTotalRecoHits(clockData, trkHits, false));
+      // T0
+      double t0 = std::numeric_limits<float>::signaling_NaN();
+      if( auto const& t0ref = fmt0pandora.at(p_pfp.key())) t0 = t0ref.ref().Time();
+
+      int goodHits=0;
+      int countE=0, countW=0;
+      
+      // These counters are used to determine if track is CC-E, EE, EW, CC-W, WE, WW
+      // depending on the track type, the Top CRT uses the appropriate position corrections
+      std::vector<double> ht;
+      std::vector<geo::Point_t> positionVector;
+      for(auto const& [trkHit, trkHitMeta]: util::zip(trkHits, trkHitMetas)){
+        bool badhit = (trkHitMeta->Index() == std::numeric_limits<unsigned int>::max()) ||
+                    (!track.HasValidPoint(trkHitMeta->Index()));
+        if(badhit) continue;
+        geo::Point_t loc = track.LocationAtPoint(trkHitMeta->Index());
+        positionVector.push_back(loc);
+        ht.push_back(trkHit->PeakTime());
+        goodHits++;
+        if(trkHit->WireID().TPC==0 || trkHit->WireID().TPC==1) countE++;
+        else countW++;
+      }
+      if(goodHits<fMinimumGoodHits) {
+        mf::LogDebug("CRTT0Tagging:")<<"Track does not have the minimal requirements of good hits: "<<goodHits;
+        continue;
+      }
+      int trackType=-1;
+      int const cryo = trkHits[0]->WireID().Cryostat;
+      if(countW!=0 && countE!=0 && cryo==0) trackType=0; //CCEast
+      else if(countW!=0 && countE==0 && cryo==0) trackType=2; //East-West
+      else if(countW==0 && countE!=0 && cryo==0) trackType=1; //East-East      
+      else if(countW!=0 && countE!=0 && cryo==1) trackType=3; //CCWest
+      else if(countW!=0 && countE==0 && cryo==1) trackType=5; //West-West
+      else if(countW==0 && countE!=0 && cryo==1) trackType=4; //West-East     
+      icarus::crt::dataTools::TopCRTCorrectionMap TopCRTCorrection;
+      if(trackType==0) TopCRTCorrection=fTopCRTTransformations.EastCC;
+      else if(trackType==1) TopCRTCorrection=fTopCRTTransformations.EE;
+      else if(trackType==2) TopCRTCorrection=fTopCRTTransformations.EW;
+      else if(trackType==3) TopCRTCorrection=fTopCRTTransformations.WestCC;
+      else if(trackType==4) TopCRTCorrection=fTopCRTTransformations.WE;
+      else if(trackType==5) TopCRTCorrection=fTopCRTTransformations.WW;
+      
+      std::vector<icarus::crt::CandCRT> crtCands;
+      for(art::Ptr<CRTHit> const& p_crthit: CRTHitList){
+        const CRTHit &crtHit = *p_crthit;
+        double crtTime=crtHit.ts1_ns/1e3;
+        // If the Track has a Pandora T0, this is also used to look for compatible CRT Hits 
+        if(!isnan(t0)){
+          if(fabs(t0-crtHit.ts1_ns)>fMaximumDeltaT) continue;
+        }
+        icarus::crt::DriftedTrack thisDriftedTrack = fMatchingAlg.DriftTrack(trkHits, trkHitMetas, fGeometryService, detProp, clockData, crtTime, track, 0);    
+        if(thisDriftedTrack.outbound>0) continue;
+        icarus::crt::PCAResults driftedPCAResults=fMatchingAlg.PCAfit(thisDriftedTrack.sp);
+        icarus::crt::TranslationVector translVector = {driftedPCAResults.eigenVector1, driftedPCAResults.mean};
+        int crtSys=fCrtUtils.MacToTypeCode(crtHit.feb_id[0]);
+        if(crtSys==2) continue; // lets discard Bottom CRT Hits for the moment
+
+        geo::Point_t delta = {std::numeric_limits<float>::signaling_NaN(), std::numeric_limits<float>::signaling_NaN(), std::numeric_limits<float>::signaling_NaN()};
+
+        double crtDistance=std::numeric_limits<float>::signaling_NaN();
+
+        icarus::crt::CRTPlane thisCRTPlane = fMatchingAlg.DeterminePlane(crtHit);
+        icarus::crt::CrossingPoint crossPoint = fMatchingAlg.DetermineProjection(translVector, thisCRTPlane);
+
+        geo::Point_t CRTHitCoordinate = {crtHit.x_pos, crtHit.y_pos, crtHit.z_pos};
+
+        if(fData && fTopCRTTransformations.imported){ // Realignment only applies to Data, not MC
+          if(crtSys==0){
+            CRTHitCoordinate = icarus::crt::dataTools::ApplyTransformation(crtHit, TopCRTCorrection, fTopCRTCenterMap);
+          }
+        }
+        delta=CRTHitCoordinate-crossPoint;
+        crtDistance=std::hypot(delta.X(), delta.Y(), delta.Z());
+        if(crtDistance>fMaximalCRTDistance) continue;
+        icarus::crt::CandCRT thisCrtCand={crtHit,p_crthit, thisCRTPlane.first, crtDistance, delta, crossPoint};
+        crtCands.push_back(thisCrtCand);
+      } // End of CRT Hit loop
+      if(crtCands.empty()) {
+        mf::LogDebug("CRTT0Tagging:")<<"No Good CRT match candidates for this track";
+        continue;
+      }
+      auto minElementIt = std::min_element(crtCands.begin(), crtCands.end(), [](const icarus::crt::CandCRT& a, const icarus::crt::CandCRT& b) {
+        return a.distance < b.distance;
+      });
+      icarus::crt::CandCRT& bestCrtCand=*minElementIt;
+      if(bestCrtCand.distance<=fGoodCandidateDistance){ 
+        int matchedSys=fCrtUtils.MacToTypeCode(bestCrtCand.CRThit.feb_id[0]);
+        if(matchedSys==2) continue; // lets discard Bottom CRT Hits for the moment
+        bool truthFound=false;
+        bool trueMatch=false;
+        bool truthIsNu=false;
+        int trueG4TrackId=std::numeric_limits<int>::lowest();
+        int truePdg=std::numeric_limits<int>::lowest();
+      
+        if(!fData && !fSkipTruth){
+          std::vector<int> crtTracks, crtPdgs;
+          std::pair<int,double> thisMatch=std::make_pair((int)bestCrtCand.CRThit.feb_id[0], bestCrtCand.CRThit.ts1_ns);
+          auto const itMatch = crtParticleMap.find(thisMatch);
+          if(itMatch != crtParticleMap.end()){
+            for(int trackID: itMatch->second){
+              crtTracks.push_back(trackID);
+              crtPdgs.push_back(particleMap.at(trackID)->PdgCode());
+            }
+          }
+          for(int trackID: crtTracks){
+            if(trackID == trueTrackId){
+              trueMatch=true;
+              break;
+            }
+          }
+          trueG4TrackId=abs(RecoUtils::TrueParticleIDFromTotalRecoHits(clockData, trkHits, false));
+          if(trueG4TrackId==0){
+            trueG4TrackId=std::numeric_limits<int>::lowest();
+          } else {
+            truthFound=true;
+            truthIsNu=isNuMap.at(trueG4TrackId);
+            truePdg=particleMap.at(trueG4TrackId)->PdgCode();
+          }
+        }
+        icarus::crt::DriftedTrack thisMatchedDriftedTrack = fMatchingAlg.DriftTrack(trkHits, trkHitMetas, fGeometryService, detProp, clockData, bestCrtCand.CRThit.ts1_ns/1e3, track, 0);
+        icarus::crt::PCAResults driftedMatchedPCAResults=fMatchingAlg.PCAfit(thisMatchedDriftedTrack.sp);
+        fEvent=e.event();
+        fRun=e.run();
+        fCrtRegion=bestCrtCand.CRThit.plane;
+        fCrtSys=matchedSys;
+        fCryo=cryo;
+        fTrackLength=track.Length();
+        fFirstEigenValue=driftedMatchedPCAResults.eigenValue1;
+        fSecondEigenValue=driftedMatchedPCAResults.eigenValue2;
+        fThirdEigenValue=driftedMatchedPCAResults.eigenValue3;
+        fTrackCrtDistance=bestCrtCand.distance;
+        fTrackCrtDeltaX=bestCrtCand.delta.X();
+        fTrackCrtDeltaY=bestCrtCand.delta.Y();
+        fTrackCrtDeltaZ=bestCrtCand.delta.Z();
+        fCrtX=bestCrtCand.CRThit.x_pos;
+        fCrtY=bestCrtCand.CRThit.y_pos;
+        fCrtZ=bestCrtCand.CRThit.z_pos;
+        fCrossPointX=bestCrtCand.crossPoint.X();
+        fCrossPointY=bestCrtCand.crossPoint.Y();
+        fCrossPointZ=bestCrtCand.crossPoint.Z();    
+        fTrueMatch=trueMatch;
+        fTree->Fill();
+
+        sbn::crt::CRTTaggingTrackFit trackFit = sbn::crt::CRTTaggingTrackFit::pca;
+        sbn::crt::CRTTaggingMethod matchMethod = sbn::crt::CRTTaggingMethod::crtHits;
+
+        mf::LogInfo("CRTT0Tagging")
+	        <<"Matched CRT time = "<<bestCrtCand.CRThit.ts1_ns/1e3<<" [us] to track "<<track.ID()<<" with projection-hit distance = "<<bestCrtCand.distance<<" Track T0 "<<t0
+	        <<"\nMatched CRT hit plane: "<<bestCrtCand.CRThit.plane<<" xpos "<<bestCrtCand.CRThit.x_pos<<" ypos "<<bestCrtCand.CRThit.y_pos<<" zpos "<<bestCrtCand.CRThit.z_pos
+          <<"\nDelta: X "<<bestCrtCand.delta.X()<<" Y "<<bestCrtCand.delta.Y()<<" Z "<<bestCrtCand.delta.Z();
+
+        t0col->push_back(anab::T0(bestCrtCand.CRThit.ts1_ns, track.ID(), matchedSys, bestCrtCand.CRThit.plane,bestCrtCand.distance));
+        art::Ptr<anab::T0> const newT0ptr = makeT0ptr(t0col->size()-1); // index of the last T0
+        trackAssn->addSingle(trkPtr, newT0ptr);
+        t0CrtHitAssn->addSingle(bestCrtCand.ptrCRThit, newT0ptr);
+        pfpAssn->addSingle(p_pfp, newT0ptr);
+
+        sbn::crt::CRTHitT0TaggingInfo matchInfo {bestCrtCand.distance, matchedSys, bestCrtCand.CRThit.plane, bestCrtCand.CRThit.ts1_ns, bestCrtCand.delta.X(), bestCrtCand.delta.Y(), bestCrtCand.delta.Z(), bestCrtCand.crossPoint.X(), bestCrtCand.crossPoint.Y(), bestCrtCand.crossPoint.Z(), bestCrtCand.plane, trackFit, matchMethod};        
+        matchInfoCol->push_back(matchInfo);
+        
+        art::Ptr<sbn::crt::CRTHitT0TaggingInfo> const newMatchInfoPtr = makeMatchInfoPtr(matchInfoCol->size()-1); // index of the last CRTHitT0TaggingInfo
+        trackMatchInfoAssn->addSingle(trkPtr, newMatchInfoPtr);
+        pfpMatchInfoAssn->addSingle(p_pfp, newMatchInfoPtr);
+
+        if(!fData && !fSkipTruth){
+          sbn::crt::CRTHitT0TaggingTruthInfo truthInfo {truthFound, trueMatch, trueG4TrackId, truePdg, truthIsNu};        
+          matchInfoTruthCol->push_back(truthInfo);
+          art::Ptr<sbn::crt::CRTHitT0TaggingTruthInfo> const newMatchTruthInfoPtr = makeMatchTruthInfoPtr(matchInfoTruthCol->size()-1); // index of the last CRTHitT0TaggingTruthInfo
+          truthAssn->addSingle(newMatchInfoPtr, newMatchTruthInfoPtr);
+        }
+
+      }
+	  } // End of Track Loop
+	} // End of Cryo Loop
+  e.put(std::move(t0col));
+  e.put(std::move(trackAssn));
+  e.put(std::move(t0CrtHitAssn));
+  e.put(std::move(matchInfoCol));
+  e.put(std::move(pfpAssn));
+  e.put(std::move(trackMatchInfoAssn));
+  e.put(std::move(pfpMatchInfoAssn));
+  e.put(std::move(matchInfoTruthCol));
+  e.put(std::move(truthAssn));
+}
+
+DEFINE_ART_MODULE(CRTT0Tagging)
diff --git a/icaruscode/CRT/CRTUtils/CRTMatchingUtils.cxx b/icaruscode/CRT/CRTUtils/CRTMatchingUtils.cxx
new file mode 100644
index 000000000..aa9dc4a30
--- /dev/null
+++ b/icaruscode/CRT/CRTUtils/CRTMatchingUtils.cxx
@@ -0,0 +1,203 @@
+/**
+ * @file   icaruscode/CRT/CRTUtils/CRTMatchingUtils.cxx
+ * @author Francesco Poppi (poppi@bo.infn.it)
+ * @date   January 2025
+ */
+
+#include "icaruscode/CRT/CRTUtils/CRTMatchingUtils.h"
+
+#include <string>
+#include <limits>
+#include <stdexcept>
+#include <sstream>
+#include <fstream>
+#include <cetlib/search_path.h>
+#include "messagefacility/MessageLogger/MessageLogger.h"
+#include "TMatrixD.h"
+#include "TMatrixDSym.h"
+#include "TMatrixDEigen.h"
+#include "TMatrixDSymEigen.h"
+#include "larcorealg/Geometry/TPCGeo.h"
+#include "larcorealg/Geometry/PlaneGeo.h"
+namespace icarus::crt{
+    
+    CRTMatchingAlg::CRTMatchingAlg(const fhicl::ParameterSet& pset)
+    {
+        this->reconfigure(pset);
+        return;
+    }
+
+    CRTMatchingAlg::CRTMatchingAlg() = default;
+
+    void CRTMatchingAlg::reconfigure(const fhicl::ParameterSet& pset)
+    {
+        fAllowedOffsetCM = pset.get<double>("AllowedOffsetCM", 1.57);
+        return;
+    }
+
+    PCAResults CRTMatchingAlg::PCAfit (std::vector<geo::Point_t> const& sp)
+    {
+        int min=0;
+        int max=sp.size();
+        int size=max-min;
+        double xavg=0, yavg=0, zavg=0;
+        for(int k=min; k<max; k++) {
+            xavg+=sp[k].X(); yavg+=sp[k].Y(); zavg+=sp[k].Z();
+        }
+        xavg=xavg/size; yavg=yavg/size; zavg=zavg/size;
+        double x2=0, y2=0, z2=0, xiy=0, xiz=0, yiz=0;
+        for(int k=min; k<max; k++){
+            double xadj = sp[k].X() - xavg;
+            double yadj = sp[k].Y() - yavg;
+            double zadj = sp[k].Z() - zavg;
+            x2+=xadj*xadj; y2+=yadj*yadj; z2+=zadj*zadj;
+            xiy+=xadj*yadj; xiz+=xadj*zadj; yiz+= yadj*zadj;
+        }//end loop calculating covariance matrix elements
+        x2=x2/size; y2=y2/size; z2=z2/size; 
+        xiy=xiy/size; xiz=xiz/size; yiz=yiz/size;
+        //Now covariance elements are ready, make covariance matrix and do PCA analysis
+        TMatrixDSym covmat(3);
+        covmat(0, 0) = x2; covmat(1, 1) = y2; covmat(2, 2) = z2;
+        covmat(0, 1) = xiy; covmat(1, 0) = xiy;
+        covmat(0, 2) = xiz; covmat(2, 0) = xiz;
+        covmat(1, 2) = yiz; covmat(2, 1) = yiz;
+        TMatrixDSymEigen thiscov(covmat);
+        TVectorD this_eval = thiscov.GetEigenValues();
+        TMatrixD this_evec = thiscov.GetEigenVectors();
+
+        if(this_eval.NonZeros()!=3) throw std::logic_error("CRTMatchingAlg::PCAfit: evals/evects wrong size!");
+        geo::Point_t mean = {xavg, yavg, zavg};        
+        geo::Vector_t first = {this_evec[0][0], this_evec[1][0], this_evec[2][0]};
+        geo::Vector_t second = {this_evec[0][1], this_evec[1][1], this_evec[2][1]};
+        geo::Vector_t third = {this_evec[0][2], this_evec[1][2], this_evec[2][2]};
+
+        PCAResults thPCAResult = {first, second, third, this_eval[0], this_eval[1], this_eval[2], mean};
+        return thPCAResult;
+    }
+
+    CRTPlane CRTMatchingAlg::DeterminePlane(sbn::crt::CRTHit const& CRThit)
+    {
+        int Plane;
+        double Pos;
+        switch (CRThit.plane) {
+          case 30:
+            Plane=0;
+            Pos=CRThit.y_pos;
+            break;
+          case 31: case 32:
+          case 40: case 41: case 42: case 43: case 44: case 45:
+            Plane=1;
+            Pos=CRThit.x_pos;
+            break;
+          default:
+            Plane=2;
+            Pos=CRThit.z_pos;
+            break;
+        } // switch
+        return std::make_pair(Plane, Pos);
+    }
+
+    CrossingPoint CRTMatchingAlg::TranslatePointTo(geo::Vector_t dirVector, geo::Point_t meanPos, CRTPlane CRTwall)
+    {
+        
+        if(dirVector.X()==0) throw std::invalid_argument("CRTMatchingAlg::TranslatePointTo: Cosine Director is null");
+        return meanPos + dirVector * (CRTwall.second - meanPos.X()) / dirVector.X();
+    }
+
+    TranslationVector CRTMatchingAlg::RotateToLocalCRTPlane(const TranslationVector& transl, CRTPlane CRTwall)
+    {
+        geo::Vector_t rotatedDir;
+        geo::Point_t rotatedMean;
+
+        switch(CRTwall.first) {
+            case 0: // Fixed Coordinate: Y
+                rotatedDir={transl.dir.Y(), transl.dir.X(), transl.dir.Z()};
+                rotatedMean={transl.mean.Y(), transl.mean.X(), transl.mean.Z()};
+                break;
+            case 1: // Fixed Coordinate: X
+                rotatedDir={transl.dir.X(), transl.dir.Y(), transl.dir.Z()};
+                rotatedMean={transl.mean.X(), transl.mean.Y(), transl.mean.Z()};
+                break;
+            case 2: // Fixed Coordinate: Z
+                rotatedDir={transl.dir.Z(), transl.dir.Y(), transl.dir.X()};
+                rotatedMean={transl.mean.Z(), transl.mean.Y(), transl.mean.X()};
+                break;
+            default:
+                throw std::invalid_argument("CRTMatchingAlg::RotateToLocalCRTPlane(): invalid plane");
+        }
+        return {rotatedDir, rotatedMean};
+    }
+
+    CrossingPoint CRTMatchingAlg::RotateFromLocalCRTPlane(CrossingPoint crossPointCRT, CRTPlane CRTWall)
+    {
+        switch(CRTWall.first) {
+            case 0: return {crossPointCRT.Y(), crossPointCRT.X(), crossPointCRT.Z()}; // Plane at Y e.g. Top CRT Horizontal Plane
+            case 1: return {crossPointCRT.X(), crossPointCRT.Y(), crossPointCRT.Z()}; // Plane at X e.g. Side CRT West, East, Top CRT
+            case 2: return {crossPointCRT.Z(), crossPointCRT.Y(), crossPointCRT.X()}; // Plane at Z e.g. Side CRT South, North, Top CRT 
+            default:
+                throw std::invalid_argument("CRTMatchingAlg::RotateFromLocalCRTPlane(): invalid plane");
+        }
+    }
+
+    CrossingPoint CRTMatchingAlg::DetermineProjection(const TranslationVector& dir, CRTPlane CRTWall)
+    {
+        TranslationVector directionPlaneCoordinate = CRTMatchingAlg::RotateToLocalCRTPlane(dir, CRTWall);
+        CrossingPoint crossingPointPlaneCoordinate = CRTMatchingAlg::TranslatePointTo(directionPlaneCoordinate.dir, directionPlaneCoordinate.mean, CRTWall);
+        return CRTMatchingAlg::RotateFromLocalCRTPlane(crossingPointPlaneCoordinate, CRTWall);
+    }
+
+    TrackBarycenter CRTMatchingAlg::GetTrackBarycenter (std::vector<geo::Point_t> sp, std::vector<double> hw)
+    {
+        double average_x_charge=0, average_y_charge=0, average_z_charge=0, total_charge=0;
+        bool isGood;
+        for (unsigned i = 0; i < sp.size(); i++) {
+            if(isnan(sp[i].Z())) continue;
+
+            if(sp[i].Z()>-1000 && sp[i].Z()<1000){
+                average_z_charge+=sp[i].Z()*hw[i];
+                average_y_charge+=sp[i].Y()*hw[i];
+                average_x_charge+=sp[i].X()*hw[i];
+                total_charge+=hw[i];
+            }
+        }
+        average_z_charge=average_z_charge/total_charge;
+        average_y_charge=average_y_charge/total_charge;
+        average_x_charge=average_x_charge/total_charge;
+        if(total_charge==0) isGood=false;
+        else isGood=true;
+        TrackBarycenter ThisTrackBary = {average_x_charge, average_y_charge, average_z_charge, isGood};
+        return ThisTrackBary;
+    }
+
+    DriftedTrack CRTMatchingAlg::DriftTrack(const std::vector<art::Ptr<recob::Hit>>& trkHits, const std::vector<const recob::TrackHitMeta*>& trkHitMetas, const geo::GeometryCore *GeometryService, detinfo::DetectorPropertiesData const& detProp, detinfo::DetectorClocksData  const& detClock, double time, const recob::Track& tpcTrack, int maxOutBoundPoints) const
+    {
+        int outBound=0;
+        std::vector<double> recI;
+        std::vector<geo::Point_t> driftedPositionVector;
+        for(size_t i=0; i<trkHits.size(); i++){
+            bool badhit = !trkHitMetas[i] || (trkHitMetas[i]->Index() == std::numeric_limits<unsigned int>::max()) ||
+                        (!tpcTrack.HasValidPoint(trkHitMetas[i]->Index()));
+            if(badhit) continue;
+            geo::Point_t loc = tpcTrack.LocationAtPoint(trkHitMetas[i]->Index());
+            const geo::TPCGeo& tpcGeo = GeometryService->TPC(trkHits[i]->WireID());
+            double vDrift=detProp.DriftVelocity();
+
+            // recoX: distance of the hit charge from the plane
+            //  * `trkHits[i]->PeakTime()-fTickAtAnode`: TPC ticks from the trigger to when charge gets to plane
+            //  * `time`: t0 [CRT or Flash] (with respect to trigger) in TPC ticks
+            //  * difference is how many ticks passed from the track to when charge gets to plane: drift ticks
+            double recoX=(trkHits[i]->PeakTime()-detClock.Time2Tick(detClock.TriggerTime())-time/detClock.TPCClock().TickPeriod())*detClock.TPCClock().TickPeriod()*vDrift;
+            double plane=tpcGeo.PlanePtr(trkHits[i]->WireID().Plane)->GetCenter().X();
+            double X=plane-tpcGeo.DriftDir().X()*recoX;
+            double AllowedRel = 1+fAllowedOffsetCM / abs(tpcGeo.GetCathodeCenter().X()-tpcGeo.PlanePtr(trkHits[i]->WireID().Plane)->GetCenter().X());
+            if(!tpcGeo.ActiveBoundingBox().ContainsX((X), AllowedRel)) outBound++;
+            if(outBound>maxOutBoundPoints) {
+                return {{},{},maxOutBoundPoints+1};
+            }
+            driftedPositionVector.push_back({X, loc.Y(), loc.Z()});
+            recI.push_back(trkHits[i]->Integral());
+        }
+        DriftedTrack thisDriftedTrack = {driftedPositionVector, recI, outBound};
+        return thisDriftedTrack;
+    }
+}
\ No newline at end of file
diff --git a/icaruscode/CRT/CRTUtils/CRTMatchingUtils.h b/icaruscode/CRT/CRTUtils/CRTMatchingUtils.h
new file mode 100644
index 000000000..fdb607d08
--- /dev/null
+++ b/icaruscode/CRT/CRTUtils/CRTMatchingUtils.h
@@ -0,0 +1,129 @@
+#ifndef ICARUSCODE_CRT_CRTUTILS_CRTMATCHINGUTILS_H
+#define ICARUSCODE_CRT_CRTUTILS_CRTMATCHINGUTILS_H
+
+/**
+ * @file   icaruscode/CRT/CRTUtils/CRTMatchingUtils.h
+ * @author Francesco Poppi (poppi@bo.infn.it)
+ * @date   January 2025
+ */
+
+#include <vector>
+#include <utility> // std::pair
+#include <map>
+
+#include "fhiclcpp/ParameterSet.h"
+#include "canvas/Persistency/Common/Ptr.h" 
+// LArSoft headers
+#include "lardataobj/RecoBase/Hit.h"
+#include "lardataobj/RecoBase/Track.h"
+#include "lardataobj/RecoBase/TrackHitMeta.h"
+#include "larcorealg/Geometry/GeometryCore.h"
+#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
+#include "lardata/DetectorInfoServices/DetectorClocksService.h"
+#include "lardataalg/DetectorInfo/DetectorPropertiesData.h"
+#include "lardataalg/DetectorInfo/DetectorClocksData.h"
+#include "lardataalg/DetectorInfo/DetectorClocks.h"
+// Data product headers
+#include "sbnobj/Common/CRT/CRTHit.hh"
+#include "sbnobj/Common/Trigger/ExtraTriggerInfo.h"
+#include "sbnobj/Common/CRT/CRTPMTMatching.hh"
+
+
+
+namespace icarus::crt{
+
+struct TrackBarycenter
+{
+    double BarX; // Track Barycenter X coordinate
+    double BarY; // Track Barycenter Y coordinate
+    double BarZ; // Track Barycenter Z coordinate
+    bool isGood; // Track Barycenter quality
+};
+
+struct DriftedTrack
+{
+    std::vector<geo::Point_t> sp; // Drifted space points
+    std::vector<double> spi; // Drifted Track Hit Points integral
+    int outbound; // Number of hit points out of the logical volume of the TPC
+    //double drifted_startx; 
+    //double drifted_endx;
+};
+
+struct TranslationVector{
+    geo::Vector_t dir;
+    geo::Point_t mean;
+};
+
+struct PCAResults{
+    geo::Vector_t eigenVector1; // First EigenVector
+    geo::Vector_t eigenVector2; // Second EigenVector
+    geo::Vector_t eigenVector3; // Third EigenVector
+    double eigenValue1; // First EigenValue
+    double eigenValue2; // Second EigenValue
+    double eigenValue3; // Third EigenValue
+    geo::Point_t mean; // Mean X,Y,Z coordinates
+};
+
+using CrossingPoint = geo::Point_t;
+
+struct CandCRT{
+    sbn::crt::CRTHit CRThit;
+    art::Ptr<sbn::crt::CRTHit> ptrCRThit;
+    int plane;
+    double distance;
+    geo::Point_t delta;
+    CrossingPoint crossPoint;
+};
+
+/// CRTPlane corresponds to a pair of:
+/// 1st CRT fixed coordinate plane (e.g. 0 is X coordinate)
+/// 2nd value of the  CRT fixed coordinate plane (e.g. for Top CRT Horizontal this value is ~618 cm.
+using CRTPlane = std::pair<int,double>;
+
+class CRTMatchingAlg {
+public:
+
+    explicit CRTMatchingAlg(const fhicl::ParameterSet& pset);
+    CRTMatchingAlg();
+
+    void reconfigure(const fhicl::ParameterSet& pset);
+    
+    /// This function runs the PCA analysis on the spatial coordinates vector and return the PCAResults.
+    static PCAResults PCAfit (std::vector<geo::Point_t> const& sp);
+    
+    /// This function determines the coordinate in which the CRT module is constant.
+    /// 0 for fixed Y, 1 for fixed X, 2 for fixed Z,
+    /// e.g. in the Top CRT Horizontal Plane, the Y coordinate is fixed, in the Side CRT West Walll, the X Coordinate is fixed, ...
+    CRTPlane DeterminePlane(sbn::crt::CRTHit const& CRThit);
+
+    /// This function evaluates the Track Crossing point onto the CRT plane considered.
+    /// dir is the track direction in the CRTWall reference system.
+    /// mean is the mean value of the track in the CRTWall reference system.
+    static CrossingPoint TranslatePointTo(geo::Vector_t dir, geo::Point_t mean, CRTPlane CRTWall);
+
+    /// This function rotate the translation vector (direction and mean position) from the TPC reference system, to the CRTWall one.
+    TranslationVector RotateToLocalCRTPlane(const TranslationVector& transl, CRTPlane CRTWall);
+
+    /// This function rotate the CrossingPoint from the CRTWall reference system, to the TPC one.
+    CrossingPoint RotateFromLocalCRTPlane(CrossingPoint crossPointCRT, CRTPlane CRTWall);
+
+    /// This function returns the "correct" predicted crossing point given the cosine directors of the track fitted with PCA
+    /// and a CRTPlane (fixed coordinate and value of the coordinate).
+    CrossingPoint DetermineProjection(const TranslationVector& dir, CRTPlane CRTWall);
+    
+    /// This function evaluates the TrackBarycenter with weighted mean. posVector is a vector of track spacepoints
+    /// and w is the weight (e.g. one could use integral or peak amplitude or more...).
+    /// The entries in the i-th index of posVector and hw vectors must correspond to the same point.
+    TrackBarycenter GetTrackBarycenter(std::vector<geo::Point_t> posVector, std::vector<double> hw);
+
+    /// Function which drifts the track x coordinates assuming a time (in the trigger reference system).
+    /// The function also returns the number of hit points which would be outside of the physical boundaries of the TPCs at the considered time.
+    DriftedTrack DriftTrack(const std::vector<art::Ptr<recob::Hit>>& trkHits, const std::vector<const recob::TrackHitMeta*>& trkHitMetas, const geo::GeometryCore *GeometryService, detinfo::DetectorPropertiesData const& detProp, detinfo::DetectorClocksData  const& detClock, double time, const recob::Track& tpcTrack, int maxOutBoundPoints) const;
+
+private:
+    double fAllowedOffsetCM;
+};
+
+}
+
+#endif // CRTMATCHINGUTILS_H
\ No newline at end of file
diff --git a/icaruscode/CRT/CRTUtils/CRTT0MatchAlg.h b/icaruscode/CRT/CRTUtils/CRTT0MatchAlg.h
index 69d30bd46..fffa3abc0 100644
--- a/icaruscode/CRT/CRTUtils/CRTT0MatchAlg.h
+++ b/icaruscode/CRT/CRTUtils/CRTT0MatchAlg.h
@@ -1,13 +1,11 @@
 #ifndef CRTT0MATCHALG_H_SEEN
 #define CRTT0MATCHALG_H_SEEN
 
-
-///////////////////////////////////////////////
-// CRTT0MatchAlg.h
-//
-// Functions for CRT t0 matching
-// T Brooks (tbrooks@fnal.gov), November 2018
-///////////////////////////////////////////////
+/**
+ * @file   icaruscode/CRT/CRTUtils/CRTT0MatchAlg.h
+ * @author Francesco Poppi (poppi@bo.infn.it)
+ * @date   January 2025
+ */
 
 // framework
 #include "art/Framework/Principal/Event.h"
@@ -45,16 +43,6 @@
 //#include "icaruscode/Geometry/GeometryWrappers/TPCGeoAlg.h"
 #include "icaruscode/CRT/CRTUtils/TPCGeoUtil.h"
 
-// c++
-#include <iostream>
-#include <stdio.h>
-#include <sstream>
-#include <vector>
-#include <map>
-#include <utility>
-#include <cmath> 
-#include <memory>
-
 // ROOT
 #include "TVector3.h"
 #include "TGeoManager.h"
diff --git a/icaruscode/CRT/CRTUtils/RecoUtils.cc b/icaruscode/CRT/CRTUtils/RecoUtils.cc
index 00ca4c99e..5b4e4e7b8 100644
--- a/icaruscode/CRT/CRTUtils/RecoUtils.cc
+++ b/icaruscode/CRT/CRTUtils/RecoUtils.cc
@@ -1,6 +1,276 @@
+#include <fstream>
+#include <cetlib/search_path.h>
 #include "RecoUtils.h"
 
 
+namespace icarus::crt::dataTools{
+
+TopCRTCentersMap LoadTopCRTCenters()
+{
+  TopCRTCentersMap TopCRTCenters;
+  //  The follwing numbers have been extracted from the Top CRT modules geometry.
+  //  The numbers are respectively moduleID and X,Y,Z coordinates of the
+  //  module center (X markes the spot in the sketch below).
+  //  The CRT modules are perfect square with 8 bars per side.
+  //  The fixed coordinate (e.g. Y for Top CRT Horizontal) is returned from geometry.
+  //  The transverce coordinate is returned from the average position of P1 and P2
+  //  the average position of P1 and P3.
+  //  The transformed CRT Hits are in cm.
+  //
+  //          ---------------------------------------------------------
+  //          |      |      |      |      |      |      |      |      |       
+  //          |      |      |      |      |      |      |      |      |       
+  //          |      |      |      |      |      |      |      |      |       
+  //          ---------------------------------------------------------
+  //          |      |      |      |      |      |      |      |      |       
+  //          |      |      |      |      |      |      |      |      |       
+  //          |      |      |      |      |      |      |      |      |       
+  //          ---------------------------------------------------------
+  //          |      |      |      |      |      |      |      |      |       
+  //          |      |      |      |      |      |      |      |      |       
+  //          |      |      |      |      |      |      |      |      |       
+  //          ---------------------------------------------------------
+  //          |      |      |      |      |      |      |      |      |       
+  //          |      |      |      |  P1  |  P2  |      |      |      |       
+  //          |      |      |      |      |      |      |      |      |       
+  //          --------------------------- X ---------------------------
+  //          |      |      |      |      |      |      |      |      |       
+  //          |      |      |      |  P3  |      |      |      |      |       
+  //          |      |      |      |      |      |      |      |      |       
+  //          ---------------------------------------------------------
+  //          |      |      |      |      |      |      |      |      |       
+  //          |      |      |      |      |      |      |      |      |       
+  //          |      |      |      |      |      |      |      |      |       
+  //          ---------------------------------------------------------
+  //          |      |      |      |      |      |      |      |      |       
+  //          |      |      |      |      |      |      |      |      |       
+  //          |      |      |      |      |      |      |      |      |       
+  //          ---------------------------------------------------------
+  //          |      |      |      |      |      |      |      |      |       
+  //          |      |      |      |      |      |      |      |      |       
+  //          |      |      |      |      |      |      |      |      |       
+  //          ---------------------------------------------------------
+  //
+  TopCRTCenters = {{108, { -460.975, 617.388,-1050.61}},
+                  {109, { -460.975,617.388,-866.215}},
+                  {110, { -460.975,617.388,-681.825}},
+                  {111, { -460.975,617.388,-497.435}},
+                  {112, { -460.975,617.388,-313.045}},
+                  {113, { -460.975,617.388,-128.655}},
+                  {114, { -460.975,617.388,55.735}},
+                  {115, { -460.975,617.388,240.125}},
+                  {116, { -460.975,617.388,424.515}},
+                  {117, { -460.975,617.388,608.905}},
+                  {118, { -460.975,617.388,793.295}},
+                  {119, { -460.975,617.388,977.685}},
+                  {120, { -460.975,617.388,1162.07}},
+                  {121, { -460.975,617.388,1346.46}},
+                  {122, { -276.585,617.388,-1050.61}},
+                  {123, { -276.585,617.388,-866.215}},
+                  {124, { -276.585,617.388,-681.825}},
+                  {125, { -276.585,617.388,-497.435}},
+                  {126, { -276.585,617.388,-313.045}},
+                  {127, { -276.585,617.388,-128.655}},
+                  {128, { -276.585,617.388, 55.735}},
+                  {129, { -276.585,617.388, 240.125}},
+                  {130, { -276.585,617.388, 424.515}},
+                  {131, { -276.585,617.388, 608.905}},
+                  {132, { -276.585,617.388, 793.295}},
+                  {133, { -276.585,617.388, 977.685}},
+                  {134, { -276.585,617.388, 1162.07}},
+                  {135, { -276.585,617.388, 1346.46}},
+                  {136, { -92.195 ,617.388,-1050.61}},
+                  {137, { -92.195 ,617.388,-866.215}},
+                  {138, { -92.195 ,617.388,-681.825}},
+                  {139, { -92.195 ,617.388,-497.435}},
+                  {140, { -92.195 ,617.388,-313.045}},
+                  {141, { -92.195 ,617.388,-128.655}},
+                  {142, { -92.195 ,617.388, 55.735}},
+                  {143, { -92.195 ,617.388, 240.125}},
+                  {144, { -92.195 ,617.388, 424.515}},
+                  {145, { -92.195 ,617.388, 608.905}},
+                  {146, { -92.195 ,617.388, 793.295}},
+                  {147, { -92.195 ,617.388, 977.685}},
+                  {148, { -92.195 ,617.388, 1162.07}},
+                  {149, { -92.195 ,617.388, 1346.46}},
+                  {150, { 92.195 , 617.388, -1050.61}},
+                  {151, { 92.195 , 617.388, -866.215}},
+                  {152, { 92.195 , 617.388, -681.825}},
+                  {153, { 92.195 , 617.388, -497.435}},
+                  {154, { 92.195 , 617.388, -313.045}},
+                  {155, { 92.195 , 617.388, -128.655}},
+                  {156, { 92.195 ,617.388, 55.735}},
+                  {157, { 92.195 ,617.388, 240.125}},
+                  {158, { 92.195 ,617.388, 424.515}},
+                  {159, { 92.195 ,617.388, 608.905}},
+                  {160, { 92.195 ,617.388, 793.295}},
+                  {161, { 92.195 ,617.388, 977.685}},
+                  {162, { 92.195 ,617.388, 1162.07}},
+                  {163, { 92.195 ,617.388, 1346.46}},
+                  {164, { 276.585,617.388, -1050.61}},
+                  {165, { 276.585,617.388, -866.215}},
+                  {166, { 276.585,617.388, -681.825}},
+                  {167, { 276.585,617.388, -497.435}},
+                  {168, { 276.585,617.388, -313.045}},
+                  {169, { 276.585,617.388, -128.655}},
+                  {170, { 276.585,617.388, 55.735}},
+                  {171, { 276.585,617.388, 240.125}},
+                  {172, { 276.585,617.388, 424.515}},
+                  {173, { 276.585,617.388, 608.905}},
+                  {174, { 276.585,617.388, 793.295}},
+                  {175, { 276.585,617.388, 977.685}},
+                  {176, { 276.585,617.388, 1162.07}},
+                  {177, { 276.585,617.388, 1346.46}},
+                  {178, { 460.975,617.388, -1050.61}},
+                  {179, { 460.975,617.388, -866.215}},
+                  {180, { 460.975,617.388, -681.825}},
+                  {181, { 460.975,617.388, -497.435}},
+                  {182, { 460.975,617.388, -313.045}},
+                  {183, { 460.975,617.388, -128.655}},
+                  {184, { 460.975,617.388, 55.735}},
+                  {185, { 460.975,617.388, 240.125}},
+                  {186, { 460.975,617.388, 424.515}},
+                  {187, { 460.975,617.388, 608.905}},
+                  {188, { 460.975,617.388, 793.295}},
+                  {189, { 460.975,617.388, 977.685}},
+                  {190, { 460.975,617.388, 1162.07}},
+                  {191, { 460.975,617.388, 1346.46}},
+                  {192, { 555.265, 496.038, -1050.61}},
+                  {193, { 555.265, 496.038, -866.215}},
+                  {194, { 555.265, 496.038, -681.825}},
+                  {195, { 555.265, 496.038, -497.435}},
+                  {196, { 555.265, 496.038, -313.045}},
+                  {197, { 555.265, 496.038, -128.655}},
+                  {198, { 555.265, 496.038, 55.735}},
+                  {199, { 555.265, 496.038, 240.125}},
+                  {200, { 555.265, 496.038, 424.515}},
+                  {201, { 555.265, 496.038, 608.905}},
+                  {202, { 555.265, 496.038, 793.295}},
+                  {203, { 555.265, 496.038, 977.685}},
+                  {204, { 555.265, 496.038, 1162.07}},
+                  {205, { 555.265, 496.038, 1346.46}},
+                  {206, { -555.265, 496.038, -1050.61}},
+                  {207, { -555.265, 496.038, -866.215}},
+                  {208, { -555.265, 496.038, -681.825}},
+                  {209, { -555.265, 496.038, -497.435}},
+                  {210, { -555.265, 496.038, -313.045}},
+                  {211, { -555.265, 496.038, -128.655}},
+                  {212, { -555.265, 496.038, 55.735}},
+                  {213, { -555.265, 496.038, 240.125}},
+                  {214, { -555.265, 496.038, 424.515}},
+                  {215, { -555.265, 496.038, 608.905}},
+                  {216, { -555.265, 496.038, 793.295}},
+                  {217, { -555.265, 496.038, 977.685}},
+                  {218, { -555.265, 496.038, 1162.07}},
+                  {219, { -555.265, 496.038, 1346.46}},
+                  {220, { -460.975, 496.038, -1143.4}},
+                  {221, { -276.585, 496.038, -1143.4}},
+                  {222, { -92.195, 496.038, -1143.4}},
+                  {223, { 92.195, 496.038, -1143.4}},
+                  {224, { 276.585, 496.038, -1143.4}},
+                  // This module does not exist in reality, but exists in simulation
+                  {225, { 0, 0, 0}},
+                  {226, { -460.975, 525.038, 1533.608}},
+                  {227, { -276.585, 525.038, 1533.608}},
+                  {228, { -92.195, 525.038, 1533.608}},
+                  {229, { 92.195, 525.038, 1533.608}},
+                  {230, { 276.585, 525.038, 1533.608}},
+                  {231, { 460.975, 525.038, 1533.608}}};
+  return TopCRTCenters;
+}
+
+TransformedCRTHit AffineTransformation(double DX, double DZ,AffineTrans affine)
+{
+  double CRTX=affine.B1+DZ*affine.A12+DX*affine.A11;
+  double CRTZ=affine.B2+DZ*affine.A22+DX*affine.A21;
+  return std::make_pair(CRTX, CRTZ);
+}
+
+TopCRTTransformations LoadTopCRTTransformations()
+{
+  std::string fullFileName;
+  cet::search_path searchPath("FW_SEARCH_PATH");
+  bool imported = false;
+  searchPath.find_file("TopCrtCorrectionByTPC.txt", fullFileName);
+  if (!searchPath.find_file("TopCrtCorrectionByTPC.txt", fullFileName)) {
+      mf::LogError("CRTMatchingUtils_LoadTopCrtTransformation")
+      << "Top CRT Correction transformation file not found in FW_SEARCH_PATH.";
+  }
+  std::ifstream corrFile(fullFileName, std::ios::in);
+  std::map<int, AffineTrans> Type0Corr;
+  std::map<int, AffineTrans> Type1Corr;
+  std::map<int, AffineTrans> Type2Corr;
+  std::map<int, AffineTrans> Type3Corr;
+  std::map<int, AffineTrans> Type4Corr;
+  std::map<int, AffineTrans> Type5Corr;
+  if (!corrFile.is_open()) {
+      mf::LogError("CRTRecoUtils_LoadTopCRTTransformation")
+      << "Failed to open Top CRT Correction transformation file: " << fullFileName;
+  }
+  std::string line;
+  while (std::getline(corrFile, line)) {
+      std::istringstream iss(line);
+      std::vector<double> variables;
+      double var;
+      while (iss >> var) {
+          variables.push_back(var);
+      }
+      int pos=0;
+      AffineTrans Type0 = {variables.at(1+pos),variables.at(2+pos),variables.at(3+pos),variables.at(4+pos),variables.at(5+pos),variables.at(6+pos), variables.at(7+pos),variables.at(8+pos)};
+      pos=pos+8;
+      AffineTrans Type1 = {variables.at(1+pos),variables.at(2+pos),variables.at(3+pos),variables.at(4+pos),variables.at(5+pos),variables.at(6+pos), variables.at(7+pos),variables.at(8+pos)};
+      pos=pos+8;
+      AffineTrans Type2 = {variables.at(1+pos),variables.at(2+pos),variables.at(3+pos),variables.at(4+pos),variables.at(5+pos),variables.at(6+pos), variables.at(7+pos),variables.at(8+pos)};
+      pos=pos+8;
+      AffineTrans Type3 = {variables.at(1+pos),variables.at(2+pos),variables.at(3+pos),variables.at(4+pos),variables.at(5+pos),variables.at(6+pos), variables.at(7+pos),variables.at(8+pos)};
+      pos=pos+8;
+      AffineTrans Type4 = {variables.at(1+pos),variables.at(2+pos),variables.at(3+pos),variables.at(4+pos),variables.at(5+pos),variables.at(6+pos), variables.at(7+pos),variables.at(8+pos)};
+      pos=pos+8;
+      AffineTrans Type5 = {variables.at(1+pos),variables.at(2+pos),variables.at(3+pos),variables.at(4+pos),variables.at(5+pos),variables.at(6+pos), variables.at(7+pos),variables.at(8+pos)};
+      Type0Corr.insert({variables.at(0), Type0});
+      Type1Corr.insert({variables.at(0), Type1});
+      Type2Corr.insert({variables.at(0), Type2});
+      Type3Corr.insert({variables.at(0), Type3});
+      Type4Corr.insert({variables.at(0), Type4});
+      Type5Corr.insert({variables.at(0), Type5});
+      imported = true;
+  }
+  TopCRTTransformations LoadedTransformations={Type2Corr, Type1Corr, Type0Corr, Type4Corr, Type5Corr, Type3Corr, imported};
+  return LoadedTransformations;
+}
+
+geo::Point_t ApplyTransformation(const sbn::crt::CRTHit& crthit, const TopCRTCorrectionMap& TopCRTCorrections, const TopCRTCentersMap& TopCRTCenters){
+  double centerDX=crthit.x_pos - TopCRTCenters.at((int)crthit.feb_id[0]).X();
+  double centerDY=crthit.y_pos - TopCRTCenters.at((int)crthit.feb_id[0]).Y();
+  double centerDZ=crthit.z_pos - TopCRTCenters.at((int)crthit.feb_id[0]).Z();
+  AffineTrans thisAffine=TopCRTCorrections.at((int)crthit.feb_id[0]);
+  std::pair<double,double> transCrt;
+  double crtX=crthit.x_pos, crtY=crthit.y_pos, crtZ=crthit.z_pos;
+  switch (crthit.plane){
+    case 30:
+      transCrt=icarus::crt::dataTools::AffineTransformation(centerDX, centerDZ, thisAffine);
+      crtX=transCrt.first;
+      crtZ=transCrt.second;
+      break;
+    case 31: case 32:
+      transCrt=icarus::crt::dataTools::AffineTransformation(centerDY, centerDZ, thisAffine);
+      crtY=transCrt.first;
+      crtZ=transCrt.second; 
+      break;
+    case 33: case 34:
+      transCrt=icarus::crt::dataTools::AffineTransformation(centerDX, centerDY, thisAffine);
+      crtX=transCrt.first;
+      crtY=transCrt.second;
+      break;
+    default:
+      throw std::logic_error("TopCRTAffineTransformationError: CRT plane/region not valid!");
+      break;
+  }
+  return {crtX, crtY, crtZ};
+}
+
+}
+
 int RecoUtils::TrueParticleID(detinfo::DetectorClocksData const& clockData,
                               const art::Ptr<recob::Hit> hit, bool rollup_unsaved_ids) {
   std::map<int,double> id_to_energy_map;
@@ -25,7 +295,35 @@ int RecoUtils::TrueParticleID(detinfo::DetectorClocksData const& clockData,
   return likely_track_id;
 }
 
+std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> RecoUtils::PrepSimChannels(const std::vector<art::Ptr<sim::SimChannel>> &simchannels, const geo::GeometryCore &geom) {
+  std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> ret;
 
+  for (const art::Ptr<sim::SimChannel>& sc : simchannels) {
+    // Lookup the wire of this channel
+    raw::ChannelID_t channel = sc->Channel();
+    std::vector<geo::WireID> maybewire = geom.ChannelToWire(channel);
+    geo::WireID thisWire; // Default constructor makes invalid wire
+    if (!maybewire.empty()) thisWire = maybewire[0];
+    for (const auto &item : sc->TDCIDEMap()) {
+      for (const sim::IDE &ide: item.second) {
+        // indexing initializes empty vector
+        ret[abs(ide.trackID)].push_back({thisWire, &ide});
+      }
+    }
+  }
+  return ret;
+}
+
+std::map<int, std::vector<art::Ptr<recob::Hit>>> RecoUtils::buildTrackIDtoHitsMap(const std::vector<art::Ptr<recob::Hit>> &allHits, 
+  const detinfo::DetectorClocksData &clockData, const cheat::BackTrackerService &backtracker) {
+  std::map<int, std::vector<art::Ptr<recob::Hit>>> ret;
+  for (const art::Ptr<recob::Hit>& h: allHits) {
+    for (int ID: backtracker.HitToTrackIds(clockData, *h)) {
+      ret[abs(ID)].push_back(h);
+    }
+  }
+  return ret;
+}
 
 int RecoUtils::TrueParticleIDFromTotalTrueEnergy(detinfo::DetectorClocksData const& clockData, const std::vector<art::Ptr<recob::Hit> >& hits, bool rollup_unsaved_ids) {
   art::ServiceHandle<cheat::BackTrackerService> bt_serv;
@@ -175,4 +473,4 @@ double RecoUtils::CalculateTrackLength(const art::Ptr<recob::Track> track){
     length+=(next_point-this_point).Mag();
   }
   return length;
-}
+}
\ No newline at end of file
diff --git a/icaruscode/CRT/CRTUtils/RecoUtils.h b/icaruscode/CRT/CRTUtils/RecoUtils.h
index 3f90f4d0e..8d71e665c 100644
--- a/icaruscode/CRT/CRTUtils/RecoUtils.h
+++ b/icaruscode/CRT/CRTUtils/RecoUtils.h
@@ -32,6 +32,7 @@
 //#include "lardataobj/AnalysisBase/ParticleID.h"
 #include "larsim/MCCheater/BackTrackerService.h"
 #include "larcore/Geometry/Geometry.h"
+#include "sbnobj/Common/CRT/CRTHit.hh"
 namespace detinfo { class DetectorClocksData; }
 
 
@@ -43,39 +44,93 @@ namespace detinfo { class DetectorClocksData; }
 #include "TTree.h"
 
 
+namespace icarus::crt::dataTools{
+
+using ModuleCenter = geo::Point_t;
+
+struct AffineTrans
+{
+  double A11;
+  double A12;
+  double A21;
+  double A22;
+  double B1;
+  double B2;
+  double Accuracy;
+  double N;
+};
+
+using FebIndex_t = int;
+
+using TopCRTCorrectionMap = std::map<FebIndex_t, AffineTrans>;
+
+struct TopCRTTransformations
+{
+  TopCRTCorrectionMap EE;
+  TopCRTCorrectionMap EW;
+  TopCRTCorrectionMap EastCC;
+  TopCRTCorrectionMap WE;
+  TopCRTCorrectionMap WW;
+  TopCRTCorrectionMap WestCC;
+  bool imported;
+};
+
+using TopCRTCentersMap = std::map<FebIndex_t, ModuleCenter>;
+
+/// The transformed CRT Hits are in cm
+using TransformedCRTHit = std::pair<double, double>;
+
+/// This function loads the Top CRT modules centers.
+TopCRTCentersMap LoadTopCRTCenters();
+
+/// This function performs the affine transformation of the CRT hit points.
+/// The AffineTransformation requires input variables in cm.
+TransformedCRTHit AffineTransformation(double DX, double DZ, AffineTrans affine);
+
+/// This functions loads the Affine Transformation TXT files.
+TopCRTTransformations LoadTopCRTTransformations();
+
+/// This function transforms a CRTHit with AffineTransformationFunctions.
+geo::Point_t ApplyTransformation(const sbn::crt::CRTHit& crthit, const TopCRTCorrectionMap& TopCRTCorrections, const TopCRTCentersMap& TopCRTCenters);
+
+}
+
 namespace RecoUtils{
 
-  // Returns the geant4 ID which contributes the most to a single reco hit.  
-  // The matching method looks for true particle which deposits the most true energy in the reco hit.  
-  // If rollup_unsaved_ids is set to true, any unsaved daughter than 
-  // contributed energy to the hit has its energy included in its closest ancestor that was saved.
+// Returns the geant4 ID which contributes the most to a single reco hit.  
+// The matching method looks for true particle which deposits the most true energy in the reco hit.  
+// If rollup_unsaved_ids is set to true, any unsaved daughter than 
+// contributed energy to the hit has its energy included in its closest ancestor that was saved.
 
-  int TrueParticleID(detinfo::DetectorClocksData const& clockData, const art::Ptr<recob::Hit> hit, bool rollup_unsaved_ids=1); 
+int TrueParticleID(detinfo::DetectorClocksData const& clockData, const art::Ptr<recob::Hit> hit, bool rollup_unsaved_ids=1); 
   
-  // Returns the geant4 ID which contributes the most to the vector of hits.  
-  // The matching method looks for which true particle deposits the most true energy in the reco hits
+// Returns the geant4 ID which contributes the most to the vector of hits.  
+// The matching method looks for which true particle deposits the most true energy in the reco hits
 
-  int TrueParticleIDFromTotalTrueEnergy(detinfo::DetectorClocksData const& clockData, const std::vector<art::Ptr<recob::Hit> >& hits, bool rollup_unsaved_ids=1); 
+int TrueParticleIDFromTotalTrueEnergy(detinfo::DetectorClocksData const& clockData, const std::vector<art::Ptr<recob::Hit> >& hits, bool rollup_unsaved_ids=1); 
 
-  // Returns the geant4 ID which contributes the most to the vector of hits.  
-  // The matching method looks for which true particle contributes the most reconstructed charge to the hit selection 
-  // (the reco charge of each hit is correlated with each maximally contributing true particle and summed)
+// Returns the geant4 ID which contributes the most to the vector of hits.  
+// The matching method looks for which true particle contributes the most reconstructed charge to the hit selection 
+// (the reco charge of each hit is correlated with each maximally contributing true particle and summed)
   
-  int TrueParticleIDFromTotalRecoCharge(detinfo::DetectorClocksData const& clockData, const std::vector<art::Ptr<recob::Hit> >& hits, bool rollup_unsaved_ids=1); 
+int TrueParticleIDFromTotalRecoCharge(detinfo::DetectorClocksData const& clockData, const std::vector<art::Ptr<recob::Hit> >& hits, bool rollup_unsaved_ids=1); 
 
-  // Returns the geant4 ID which contributes the most to the vector of hits.  
-  // The matching method looks for which true particle maximally contributes to the most reco hits
+// Returns the geant4 ID which contributes the most to the vector of hits.  
+// The matching method looks for which true particle maximally contributes to the most reco hits
   
-  int TrueParticleIDFromTotalRecoHits(detinfo::DetectorClocksData const& clockData, const std::vector<art::Ptr<recob::Hit> >& hits, bool rollup_unsaved_ids=1); 
+int TrueParticleIDFromTotalRecoHits(detinfo::DetectorClocksData const& clockData, const std::vector<art::Ptr<recob::Hit> >& hits, bool rollup_unsaved_ids=1); 
 
-  // Checks if a position is within any of the TPCs in the geometry (user can define some distance buffer from the TPC walls)
+// Checks if a position is within any of the TPCs in the geometry (user can define some distance buffer from the TPC walls)
   
-  bool IsInsideTPC(TVector3 position, double distance_buffer); 
+bool IsInsideTPC(TVector3 position, double distance_buffer); 
 
-  // Calculates the total length of a recob::track by summing up the distances between adjacent traj. points
+// Calculates the total length of a recob::track by summing up the distances between adjacent traj. points
 
-  double CalculateTrackLength(const art::Ptr<recob::Track> track); 
+double CalculateTrackLength(const art::Ptr<recob::Track> track); 
   
+std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> PrepSimChannels(const std::vector<art::Ptr<sim::SimChannel>> &simchannels, const geo::GeometryCore &geo);
+
+std::map<int, std::vector<art::Ptr<recob::Hit>>> buildTrackIDtoHitsMap(const std::vector<art::Ptr<recob::Hit>> &allHits, const detinfo::DetectorClocksData &clockData, const cheat::BackTrackerService &backtracker);
 }
 
 #endif
diff --git a/icaruscode/CRT/CRTUtils/crtmatchingalg_icarus.fcl b/icaruscode/CRT/CRTUtils/crtmatchingalg_icarus.fcl
new file mode 100644
index 000000000..94e207321
--- /dev/null
+++ b/icaruscode/CRT/CRTUtils/crtmatchingalg_icarus.fcl
@@ -0,0 +1,33 @@
+BEGIN_PROLOG
+
+icarus_crtmatchingalg:
+{
+    AllowedOffsetCM:        1.57    # When trying to determine if a track and a time are compatible, you reconstructruct the track at that time
+                                    # you then check if the reconstructed track is out of the boundaries of the cathode and of the anode,
+                                    # this variable gives you the maximum excess of the track wrt cathode or anode to consider the track/time pair.
+                                    # Since we are allowing for MaximumDeltaT of 10 us (see above), this corresponds to a maximum excess
+                                    # of 10 us * 0.157 cm/us = 1.57 cm.
+                                    # As this is phrased, one could think of this as the two variables are related, that is not necessarily
+                                    # the case, but is still a good guess.
+}
+
+
+icarus_crtt0tagging_base:
+{
+    module_type:         "icaruscode/CRT/CRTT0Tagging"
+    CrtHitModuleLabel:  "crthit"
+    TPCTrackLabel:       ["pandoraTrackGausCryoE", "pandoraTrackGausCryoW"]     # Track producer module label
+    TRKHMLabel:  ["pandoraTrackGausCryoE", "pandoraTrackGausCryoW"]     # TrackHit Metadata producer module label
+    PFParticleLabel: 	   ["pandoraGausCryoE", "pandoraGausCryoW"] # PFParticle producer module label
+    HitLabel:   ["cluster3DCryoE" , "cluster3DCryoW"]
+    MatchingAlg:    @local::icarus_crtmatchingalg
+    MinimalTrackLength:     40.     # Minimal TrackLength to match with a CRT
+    MinimumGoodHits:        5       # Minimum number of good hits to perform the fit
+    MaximalCRTDistance:     300.    # Maximal distance between CRT Hit Candidate and Track Projection
+    GoodCandidateDistance:  100.    # 96 cm maximizes EfficiencyXPurity (both > 92%) 
+    MaximumDeltaT:          10000.     # Maximal Time difference between a T0 tagged track and CRT Hit time for the combination to be considered acceptable. [ns]
+
+}
+
+
+END_PROLOG
diff --git a/icaruscode/CRT/crtt0tagging.fcl b/icaruscode/CRT/crtt0tagging.fcl
new file mode 100644
index 000000000..c3dfd551a
--- /dev/null
+++ b/icaruscode/CRT/crtt0tagging.fcl
@@ -0,0 +1,38 @@
+#include "crtmatchingalg_icarus.fcl"
+
+BEGIN_PROLOG
+
+icarus_crtt0tagging_data:
+{
+   @table::icarus_crtt0tagging_base
+   isData:                 true    # Data
+   skipTruth:              false   # skipTruth
+
+   module_type:         "CRTT0Tagging"
+   CrtHitModuleLabel:  "crthit"
+   TPCTrackLabel:       ["pandoraTrackGausCryoE", "pandoraTrackGausCryoW"]     # Track producer module label
+   PFParticleLabel: 	   ["pandoraGausCryoE", "pandoraGausCryoW"] # PFParticle producer module label
+   TRKHMLabel:  ["pandoraTrackGausCryoE", "pandoraTrackGausCryoW"]     # TrackHit Metadata producer module label
+   HitLabel:   ["cluster3DCryoE" , "cluster3DCryoW"]
+   MatchingAlg: @local::icarus_crtmatchingalg
+}
+
+icarus_crtt0taggingmc:
+{
+   @table::icarus_crtt0tagging_base
+   isData:                 false   # MC
+   skipTruth:              false   # skipTruth
+
+   module_type:         "CRTT0Tagging"
+   CrtHitModuleLabel:  "crthit"
+   TPCTrackLabel:       ["pandoraTrackGausCryoE", "pandoraTrackGausCryoW"]     # Track producer module label
+   PFParticleLabel: 	   ["pandoraGausCryoE", "pandoraGausCryoW"] # PFParticle producer module label
+   TRKHMLabel:  ["pandoraTrackGausCryoE", "pandoraTrackGausCryoW"]     # TrackHit Metadata producer module label
+   HitLabel:   ["cluster3DCryoE" , "cluster3DCryoW"]
+   MatchingAlg: @local::icarus_crtmatchingalg
+   SimulationLabel: "largeant"
+   AuxDetSimProducerLabel: "genericcrt"
+   SimChannelProducer:  "daq:simpleSC"
+}
+
+END_PROLOG
diff --git a/icaruscode/CRT/crtt0tagging_icarus.fcl b/icaruscode/CRT/crtt0tagging_icarus.fcl
new file mode 100644
index 000000000..8b5399698
--- /dev/null
+++ b/icaruscode/CRT/crtt0tagging_icarus.fcl
@@ -0,0 +1,28 @@
+#include "services_common_icarus.fcl"
+#include "crtt0tagging.fcl"
+#include "rootoutput_icarus.fcl"
+
+process_name: CRTT0Tagging
+
+services:
+{  
+  message:                @local::icarus_message_services_prod_debug
+                          @table::icarus_common_services
+  TFileService:           { fileName: "crtt0tagging_hist.root" }
+} # services
+
+outputs.out1: @local::icarus_rootoutput
+
+# The 'physics' section defines and configures some modules to do work on each event.
+physics:
+{
+  producers:
+  {
+    CRTT0Tagging: @local::icarus_crtt0tagging_data
+  }
+
+  # Schedule job step(s) for execution by defining the analysis module for this job.
+  reco: [ CRTT0Tagging ]
+
+  stream1: [ out1 ]
+}
diff --git a/icaruscode/CRT/crtt0taggingmc_icarus.fcl b/icaruscode/CRT/crtt0taggingmc_icarus.fcl
new file mode 100644
index 000000000..19954ecbc
--- /dev/null
+++ b/icaruscode/CRT/crtt0taggingmc_icarus.fcl
@@ -0,0 +1,36 @@
+#include "services_common_icarus.fcl"
+#include "services_icarus_simulation.fcl"
+#include "simulationservices.fcl"
+#include "backtrackerservice.fcl"
+#include "crtt0tagging.fcl"
+#include "rootoutput_icarus.fcl"
+
+process_name: CRTT0Tagging
+
+services:
+{
+  ParticleInventoryService:   @local::standard_particleinventoryservice
+  BackTrackerService:         @local::standard_backtrackerservice       # from `backtrackerservice.fcl` (`larsim`)
+  message:                @local::icarus_message_services_prod_debug
+                          @table::icarus_common_services
+                            # Load the service that manages root files for histograms.
+  TFileService:           { fileName: "crtt0tagging_hist.root" }
+
+} # services
+services.BackTrackerService.BackTracker.SimChannelModuleLabel: "daq:simpleSC"
+
+outputs.out1: @local::icarus_rootoutput
+
+# The 'physics' section defines and configures some modules to do work on each event.
+physics:
+{
+  producers:
+  {
+    CRTT0Tagging: @local::icarus_crtt0taggingmc
+  }
+
+  # Schedule job step(s) for execution by defining the analysis module for this job.
+  reco: [ CRTT0Tagging ]
+
+  stream1: [ out1 ]
+}
diff --git a/icaruscode/IcarusObj/classes.h b/icaruscode/IcarusObj/classes.h
index 7e3b391ce..7430e079a 100644
--- a/icaruscode/IcarusObj/classes.h
+++ b/icaruscode/IcarusObj/classes.h
@@ -7,7 +7,6 @@
 #include "icaruscode/IcarusObj/OpDetWaveformMeta.h"
 #include "icaruscode/IcarusObj/PMTWaveformTimeCorrection.h"
 #include "icaruscode/IcarusObj/Hit.h"
-//#include "icaruscode/IcarusObj/CRTPMTMatching.h"
 
 #include "sbnobj/ICARUS/PMT/Trigger/Data/OpticalTriggerGate.h"
 #include "sbnobj/ICARUS/PMT/Trigger/Data/TriggerGateData.h"
@@ -16,6 +15,7 @@
 #include "lardataobj/AnalysisBase/T0.h"
 #include "lardataobj/RecoBase/OpFlash.h"
 #include "lardataobj/RecoBase/PFParticle.h"
+#include "lardataobj/RecoBase/Track.h"
 #include "lardataobj/Simulation/BeamGateInfo.h"
 #include "sbnobj/Common/CRT/CRTHit.hh"
 #include "lardataobj/RawData/OpDetWaveform.h"