diff --git a/sbndcode/CMakeLists.txt b/sbndcode/CMakeLists.txt
index 2827cd1af..f90ed87a6 100755
--- a/sbndcode/CMakeLists.txt
+++ b/sbndcode/CMakeLists.txt
@@ -50,3 +50,4 @@ add_subdirectory(SBNDCVN)
 add_subdirectory(TPCPMTBarycenterMatching)
 
 add_subdirectory(BlipRecoSBND)
+add_subdirectory(LightPropagationCorrection)
diff --git a/sbndcode/JobConfigurations/base/cafmakerjob_sbnd_data_base.fcl b/sbndcode/JobConfigurations/base/cafmakerjob_sbnd_data_base.fcl
index 9374f1e89..e9c608edd 100644
--- a/sbndcode/JobConfigurations/base/cafmakerjob_sbnd_data_base.fcl
+++ b/sbndcode/JobConfigurations/base/cafmakerjob_sbnd_data_base.fcl
@@ -156,6 +156,9 @@ physics.producers.cafmaker.TrackChi2PidLabel: "pandoraPidData"
 # Include 3D barycenter flahs matching 
 physics.producers.cafmaker.TPCPMTBarycenterMatchLabel: "tpcpmtbarycentermatching"
 
+# Include Light Propagation Correction for OpFlash timing
+physics.producers.cafmaker.CorrectedOpFlashLabel: "lightpropagationcorrection"
+
 # Blinding not needed for MC
 physics.producers.cafmaker.CreateBlindedCAF: false
 
diff --git a/sbndcode/JobConfigurations/standard/caf/cafmakerjob_sbnd_data_sce.fcl b/sbndcode/JobConfigurations/standard/caf/cafmakerjob_sbnd_data_sce.fcl
index 4ffb5ba94..0732aba99 100644
--- a/sbndcode/JobConfigurations/standard/caf/cafmakerjob_sbnd_data_sce.fcl
+++ b/sbndcode/JobConfigurations/standard/caf/cafmakerjob_sbnd_data_sce.fcl
@@ -13,6 +13,7 @@ physics.producers.cafmaker.SBNDCRTTrackMatchLabel: "crttrackmatchingSCE"
 physics.producers.cafmaker.OpT0Label: "opt0finderSCE"
 physics.producers.cafmaker.TPCPMTBarycenterMatchLabel: "tpcpmtbarycentermatchingSCE"
 physics.producers.cafmaker.CVNLabel: "cvnSCE"
+physics.producers.cafmaker.CorrectedOpFlashLabel: "lightpropagationcorrectionSCE"
 
 physics.producers.cnnid.ClusterModuleLabel:        "pandoraSCE"
 physics.producers.cnnid.PFParticleModuleLabel:     "pandoraSCE"
@@ -23,6 +24,10 @@ physics.producers.pandoraTrackRange.TrackLabel: "pandoraSCETrack"
 physics.producers.cafmaker.FlashMatchOpDetSuffixes: ["", "op", "ara", "opara"]
 physics.producers.cafmaker.FlashMatchSCECryoSuffixes: ["SCE"]
 
+
+
+
+
 physics.producers.pandoraShowerSelectionVars.PandoraLabel: "pandoraSCE"
 physics.producers.pandoraShowerSelectionVars.ShowerLabel:  "pandoraSCEShowerSBN"
 physics.producers.pandoraShowerCosmicDist.PandoraLabel: "pandoraSCE"
diff --git a/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl b/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl
index 6627a859b..3891e75cb 100644
--- a/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl
+++ b/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl
@@ -2,6 +2,7 @@
 #include "crt_calib_service.fcl"
 #include "opt0finder_sbnd_data.fcl"
 #include "sbnd_tpcpmt3dbarycentermatching_config.fcl"
+#include "sbnd_lightpropagationcorrection_config.fcl"
 #include "frameshift_sbnd_data.fcl"
 #include "standard_reco2_sbnd.fcl"
 
@@ -26,6 +27,8 @@ physics.producers:
     opt0finderSCE:  @local::sbnd_opt0_finder_data
     tpcpmtbarycentermatching: @local::TPCPMTBarycenterMatchProducer
     tpcpmtbarycentermatchingSCE: @local::TPCPMTBarycenterMatchProducerSCE
+    lightpropagationcorrection: @local::LightPropagationCorrection
+    lightpropagationcorrectionSCE: @local::LightPropagationCorrectionSCE
 
     ### shower reco for data
     pandoraShower:       @local::sbnd_incremental_pandoraModularShowerCreationData
@@ -41,7 +44,8 @@ physics.reco2: [ pandora, pandoraTrack, pandoraShower, pandoraShowerSBN, pandora
                  cvn, opt0finder, crtveto, crtspacepointmatching, crttrackmatching, tpcpmtbarycentermatching,
                  pandoraSCE, pandoraSCETrack, pandoraSCEShower, pandoraSCEShowerSBN, pandoraSCECaloData, pandoraSCEPidData,
                  cvnSCE, opt0finderSCE, tpcpmtbarycentermatchingSCE, crtspacepointmatchingSCE, crttrackmatchingSCE,
-                 caloskimCalorimetry, blipreco, frameshift]
+                 caloskimCalorimetry, blipreco, lightpropagationcorrectionSCE, frameshift]
+
 
 physics.analyzers.caloskim.G4producer: ""
 physics.analyzers.caloskim.SimChannelproducer: ""
diff --git a/sbndcode/LightPropagationCorrection/CMakeLists.txt b/sbndcode/LightPropagationCorrection/CMakeLists.txt
new file mode 100644
index 000000000..c6d74cd16
--- /dev/null
+++ b/sbndcode/LightPropagationCorrection/CMakeLists.txt
@@ -0,0 +1,38 @@
+
+set (
+  MODULE_LIBRARIES
+    art::Framework_Principal
+    art::Framework_Services_Registry
+    art_root_io::tfile_support
+    art_root_io::TFileService_service
+    art::Persistency_Common
+    art::Persistency_Provenance
+    art::Utilities
+
+    ROOT::Tree
+    ROOT::Core
+    
+    larsim::Utils
+    lardata::DetectorInfoServices_DetectorClocksServiceStandard_service
+    larcore::Geometry_Geometry_service
+    larcorealg::Geometry
+    lardataobj::RecoBase
+    lardataobj::MCBase
+    lardataobj::Simulation
+
+
+    sbndcode_OpDetReco_OpFlash_FlashFinder
+
+    sbnobj::Common_Reco
+    sbnobj::SBND_Timing
+    sbndcode_OpDetSim
+)
+
+cet_build_plugin(LightPropagationCorrection art::module SOURCE LightPropagationCorrection_module.cc LIBRARIES ${MODULE_LIBRARIES})
+cet_build_plugin(LightPropagationCorrectionAna art::module SOURCE LightPropagationCorrectionAna_module.cc LIBRARIES ${MODULE_LIBRARIES})
+
+install_fhicl()
+install_source()
+install_headers()
+add_subdirectory(job)
+
diff --git a/sbndcode/LightPropagationCorrection/LightPropagationCorrectionAna_module.cc b/sbndcode/LightPropagationCorrection/LightPropagationCorrectionAna_module.cc
new file mode 100644
index 000000000..f90b0db3b
--- /dev/null
+++ b/sbndcode/LightPropagationCorrection/LightPropagationCorrectionAna_module.cc
@@ -0,0 +1,150 @@
+////////////////////////////////////////////////////////////////////////
+// Class:       SBNDOpT0FinderAna
+// Plugin Type: analyzer (art v3_05_01)
+// File:        SBNDOpT0FinderAna_module.cc
+//
+// Generated at Mon Oct 12 13:52:37 2020 by Marco Del Tutto using cetskelgen
+// from cetlib version v3_10_00.
+////////////////////////////////////////////////////////////////////////
+
+#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/Run.h"
+#include "art/Framework/Principal/SubRun.h"
+#include "canvas/Utilities/InputTag.h"
+#include "canvas/Persistency/Common/FindManyP.h"
+#include "fhiclcpp/ParameterSet.h"
+#include "messagefacility/MessageLogger/MessageLogger.h"
+
+#include "lardataobj/RecoBase/Slice.h"
+#include "lardataobj/RecoBase/OpFlash.h"
+#include "sbnobj/Common/Reco/CorrectedOpFlashTiming.h"
+
+#include "art_root_io/TFileService.h"
+
+#include "TTree.h"
+
+#include <numeric>
+
+class LightPropagationCorrectionAna;
+
+
+class LightPropagationCorrectionAna : public art::EDAnalyzer {
+public:
+  explicit LightPropagationCorrectionAna(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.
+  LightPropagationCorrectionAna(LightPropagationCorrectionAna const&) = delete;
+  LightPropagationCorrectionAna(LightPropagationCorrectionAna&&) = delete;
+  LightPropagationCorrectionAna& operator=(LightPropagationCorrectionAna const&) = delete;
+  LightPropagationCorrectionAna& operator=(LightPropagationCorrectionAna&&) = delete;
+
+  // Required functions.
+  void analyze(art::Event const& e) override;
+  void beginJob() override;
+private:
+
+  std::string fLightPropCorrectionLabel; ///< The light propagation correction label producer
+
+  TTree* fTree; ///< The TTree for storing event information
+
+  double fOpFlashT0;
+  double fUpstreamTime_lightonly;
+  double fUpstreamTime_tpczcorr;
+  double fUpstreamTime_propcorr_tpczcorr;
+
+  unsigned int _eventID;
+  unsigned int _runID;
+  unsigned int _subrunID;
+
+};
+
+
+LightPropagationCorrectionAna::LightPropagationCorrectionAna(fhicl::ParameterSet const& p)
+  : EDAnalyzer{p}
+{
+  fLightPropCorrectionLabel = p.get<std::string>("LightPropCorrectionLabel");
+}
+
+void LightPropagationCorrectionAna::analyze(art::Event const& e)
+{
+  fOpFlashT0=-99999.;
+  fUpstreamTime_lightonly=-99999.;
+  fUpstreamTime_tpczcorr=-99999.;
+  fUpstreamTime_propcorr_tpczcorr=-99999.;
+
+  _eventID = -1;
+  _runID = -1;
+  _subrunID = -1;
+
+  std::cout << "Run: " << e.id().run() << " Sub: " << e.id().subRun() << " Evt: " << e.id().event() << std::endl;
+
+  _eventID = e.id().event();
+  _runID = e.id().run();
+  _subrunID = e.id().subRun();
+
+  // Get all the T0 objects
+  ::art::Handle<std::vector<sbn::CorrectedOpFlashTiming>> correctedopflash_h;
+  e.getByLabel(fLightPropCorrectionLabel, correctedopflash_h);
+  if(!correctedopflash_h.isValid() || correctedopflash_h->empty()) {
+    mf::LogWarning("fLightPropCorrectionLabel") << "No LightPropCorrectionLabel objects with label " << fLightPropCorrectionLabel << std::endl;
+    return;
+  }
+
+
+  std::vector<art::Ptr<sbn::CorrectedOpFlashTiming>> correctedopflash_v;
+  art::fill_ptr_vector(correctedopflash_v, correctedopflash_h);
+
+  // Get the T0->Slice association
+  art::FindManyP<recob::Slice> correctedflash_to_slices(correctedopflash_h, e, fLightPropCorrectionLabel);
+  art::FindManyP<recob::OpFlash> correctedflash_to_flashes(correctedopflash_h, e, fLightPropCorrectionLabel);
+
+  for (size_t i = 0; i < correctedopflash_v.size(); i++) {
+
+    // The T0 object
+    auto const correctedopflash =  correctedopflash_v[i];
+
+    // The associations from T0 to both slices and flashes
+    std::vector<art::Ptr<recob::Slice>> slice_v = correctedflash_to_slices.at(correctedopflash.key());
+    std::vector<art::Ptr<recob::OpFlash>> flash_v = correctedflash_to_flashes.at(correctedopflash.key());
+
+    // One T0 object is always associated to one and only one Slice and Flash
+    assert(slice_v.size() == 1);
+    assert(flash_v.size() == 1);
+    
+    std::cout << " Corrected flash time is " << correctedopflash->OpFlashT0 << std::endl;
+    std::cout << " Associated with slice id " << slice_v[0]->ID() << std::endl;
+    std::cout << "Corrected flash time light only  " << correctedopflash->UpstreamTime_lightonly << std::endl;
+    std::cout << "Corrected flash time tpc z corr " << correctedopflash->UpstreamTime_tpczcorr << std::endl;
+    std::cout << "Corrected flash time prop corr tpc z corr " << correctedopflash->UpstreamTime_propcorr_tpczcorr << std::endl;
+
+    fOpFlashT0 = correctedopflash->OpFlashT0;
+    fUpstreamTime_lightonly = correctedopflash->UpstreamTime_lightonly;
+    fUpstreamTime_tpczcorr = correctedopflash->UpstreamTime_tpczcorr;
+    fUpstreamTime_propcorr_tpczcorr = correctedopflash->UpstreamTime_propcorr_tpczcorr;
+    fTree->Fill();
+  }
+}
+
+void LightPropagationCorrectionAna::beginJob()
+{
+  // Implementation of optional member function here.
+  art::ServiceHandle<art::TFileService> tfs;
+  fTree = tfs->make<TTree>("LightPropagationCorrection", "Light Propagation Correction Tree");
+
+  fTree->Branch("eventID", &_eventID, "eventID/i");
+  fTree->Branch("runID", &_runID, "runID/i");
+  fTree->Branch("subrunID", &_subrunID, "subrunID/i");
+
+  fTree->Branch("fOpFlashT0", &fOpFlashT0, "OpFlashT0/d");
+  fTree->Branch("fUpstreamTime_lightonly", &fUpstreamTime_lightonly, "UpstreamTime_lightonly/d");
+  fTree->Branch("fUpstreamTime_tpczcorr", &fUpstreamTime_tpczcorr, "UpstreamTime_tpczcorr/d");
+  fTree->Branch("fUpstreamTime_propcorr_tpczcorr", &fUpstreamTime_propcorr_tpczcorr, "UpstreamTime_propcorr_tpczcorr/d");
+
+}
+
+DEFINE_ART_MODULE(LightPropagationCorrectionAna)
\ No newline at end of file
diff --git a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc
new file mode 100644
index 000000000..5675f8cf4
--- /dev/null
+++ b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc
@@ -0,0 +1,554 @@
+#include "LightPropagationCorrection_module.hh"
+
+namespace sbnd {
+    class LightPropagationCorrection;
+}
+
+sbnd::LightPropagationCorrection::LightPropagationCorrection(fhicl::ParameterSet const& p)
+    : EDProducer{p},
+    fReco2Label( p.get<std::string>("Reco2Label") ),
+    fOpT0FinderModuleLabel( p.get<std::string>("OpT0FinderModuleLabel") ),
+    fTPCPMTBarycenterFMModuleLabel( p.get<std::string>("TPCPMTBarycenterFMModuleLabel") ),
+    fOpFlashLabel_tpc0 ( p.get<std::string>("OpFlashLabel_tpc0") ),
+    fOpFlashLabel_tpc1 ( p.get<std::string>("OpFlashLabel_tpc1") ),
+    fSpacePointLabel( p.get<std::string>("SpacePointLabel") ),
+    fOpHitsModuleLabel( p.get<std::string>("OpHitsModuleLabel") ),
+    fSPECTDCLabel( p.get<std::string>("SPECTDCLabel") ),
+    fFlashMatchingTool( p.get<std::string>("FlashMatchingTool") ),
+    fSaveCorrectionTree( p.get<bool>("SaveCorrectionTree") ),
+    fSpeedOfLight( p.get<double>("SpeedOfLight") ),
+    fVGroupVIS( p.get<double>("VGroupVIS") ),
+    fVGroupVUV( p.get<double>("VGroupVUV") ),
+    fNuScoreThreshold( p.get<double>("NuScoreThreshold") ),
+    fFMScoreThreshold( p.get<double>("FMScoreThreshold") ),
+    fDebug( p.get<bool>("Debug", false) )
+    // 
+    // More initializers here.
+{
+    // Initialize the TimeCorrectionVector PerChannel
+    for(size_t i = 0; i < fWireReadout.NOpChannels(); ++i) {
+        fTimeCorrectionPerChannel.push_back(0.0); // Initialize with zero or any default value
+    }
+
+    for(unsigned int opch=0; opch<fWireReadout.NOpChannels(); opch++){
+        auto pdCenter = fWireReadout.OpDetGeoFromOpChannel(opch).GetCenter();
+        fOpDetID.push_back(opch);
+        fOpDetX.push_back(pdCenter.X());
+        fOpDetY.push_back(pdCenter.Y());
+        fOpDetZ.push_back(pdCenter.Z());
+    }
+
+    auto const& tpc = art::ServiceHandle<geo::Geometry>()->TPC();
+    fDriftDistance = tpc.DriftDistance();
+    fKinkDistance = 0.5*fDriftDistance*(1-fVGroupVUV/fVGroupVIS);
+    fVGroupVUV_I = 1./fVGroupVUV;
+    fVISLightPropTime = fDriftDistance/fVGroupVIS;
+
+    // Initialize flash algo for both TPCs 
+    auto const flash_algo  = p.get<std::string>("FlashFinderAlgo");
+    auto const flash_pset_tpc0 = p.get<lightana::Config_t>("AlgoConfig_tpc0");
+    auto algo_ptr_tpc0 = ::lightana::FlashAlgoFactory::get().create(flash_algo,flash_algo);
+    algo_ptr_tpc0->Configure(flash_pset_tpc0);
+    _mgr_tpc0.SetFlashAlgo(algo_ptr_tpc0);
+
+    auto const flash_pset_tpc1 = p.get<lightana::Config_t>("AlgoConfig_tpc1");
+    auto algo_ptr_tpc1 = ::lightana::FlashAlgoFactory::get().create(flash_algo,flash_algo);
+    algo_ptr_tpc1->Configure(flash_pset_tpc1);
+    _mgr_tpc1.SetFlashAlgo(algo_ptr_tpc1);
+
+    //Initialize flash geo tool
+    auto const flashgeo_pset = p.get<lightana::Config_t>("FlashGeoConfig");
+    _flashgeo = art::make_tool<lightana::FlashGeoBase>(flashgeo_pset);
+    
+    //Initialize flash t0 tool
+    auto const flasht0_pset = p.get<lightana::Config_t>("FlashT0Config");
+    _flasht0calculator = art::make_tool<lightana::FlashT0Base>(flasht0_pset);
+    
+    produces< std::vector<sbn::CorrectedOpFlashTiming> >();
+    produces<art::Assns<recob::Slice, sbn::CorrectedOpFlashTiming>>();
+    produces<art::Assns<recob::OpFlash, sbn::CorrectedOpFlashTiming>>();
+}
+
+void sbnd::LightPropagationCorrection::produce(art::Event & e)
+{
+    fEvent = e.id().event();
+    fRun = e.id().run();
+    fSubrun = e.id().subRun();
+
+    std::unique_ptr< std::vector<sbn::CorrectedOpFlashTiming> > correctedOpFlashTimes (new std::vector<sbn::CorrectedOpFlashTiming>);
+    art::PtrMaker<sbn::CorrectedOpFlashTiming> make_correctedopflashtime_ptr{e};
+    std::unique_ptr< art::Assns<recob::Slice, sbn::CorrectedOpFlashTiming>> newCorrectedOpFlashTimingSliceAssn (new art::Assns<recob::Slice, sbn::CorrectedOpFlashTiming>);
+    std::unique_ptr< art::Assns<recob::OpFlash, sbn::CorrectedOpFlashTiming>> newCorrectedOpFlashTimingOpFlashAssn (new art::Assns<recob::OpFlash, sbn::CorrectedOpFlashTiming>);
+
+    // --- Read Recob Slice
+    ::art::Handle<std::vector<recob::Slice>> sliceHandle;
+    e.getByLabel(fReco2Label, sliceHandle);
+    // Slice to OpT0Finder
+    //Get the handle for making the assns later on
+    ::art::Handle<std::vector<sbn::OpT0Finder>> opt0Handle;
+    e.getByLabel(fOpT0FinderModuleLabel, opt0Handle);
+    // Slice to TPCPMTBarycenterFM
+    ::art::Handle<std::vector<sbn::TPCPMTBarycenterMatch>> tpcpmtbarycenterfmHandle;
+    e.getByLabel(fTPCPMTBarycenterFMModuleLabel, tpcpmtbarycenterfmHandle);
+
+    //Read PFPs
+    ::art::Handle<std::vector<recob::PFParticle>> pfpHandle;
+    e.getByLabel(fReco2Label, pfpHandle);
+    //Read OpFlash Handle
+    art::Handle< std::vector<recob::OpFlash> > opflashListHandle_tpc0;
+    e.getByLabel(fOpFlashLabel_tpc0, opflashListHandle_tpc0);
+    art::Handle< std::vector<recob::OpFlash> > opflashListHandle_tpc1;
+    e.getByLabel(fOpFlashLabel_tpc1, opflashListHandle_tpc1);
+
+    art::FindManyP<sbn::OpT0Finder> slice_opt0finder_assns(sliceHandle, e, fOpT0FinderModuleLabel);
+    // Slice to TPCPMTBarycenterFM
+    art::FindManyP<sbn::TPCPMTBarycenterMatch> slice_tpcpmtbarycentermatching_assns(sliceHandle, e, fTPCPMTBarycenterFMModuleLabel);
+    // Slice to hits
+    art::FindManyP<recob::Hit> slice_hit_assns (sliceHandle, e, fReco2Label);
+    //Slice to PFParticles association
+    art::FindManyP<recob::PFParticle> slice_pfp_assns (sliceHandle, e, fReco2Label);
+    //PFP to vertex
+    art::FindManyP<recob::Vertex> pfp_vertex_assns(pfpHandle, e, fReco2Label);
+    //PFP to space points
+    art::FindManyP<recob::SpacePoint> pfp_sp_assns(pfpHandle, e, fSpacePointLabel);
+    //OpFlash to OpHit
+    art::FindManyP<recob::OpHit> flashToOpHitAssns_tpc0(opflashListHandle_tpc0, e, fOpFlashLabel_tpc0);
+    art::FindManyP<recob::OpHit> flashToOpHitAssns_tpc1(opflashListHandle_tpc1, e, fOpFlashLabel_tpc1);
+    // PFP Metadata
+    art::FindManyP<larpandoraobj::PFParticleMetadata> pfp_to_metadata(pfpHandle, e, fReco2Label);
+
+    // --- Store candidate slices
+    std::vector< art::Ptr<recob::Slice> > sliceVect;
+    art::fill_ptr_vector(sliceVect, sliceHandle);
+
+    //Vector for recob PFParticles
+    std::vector<art::Ptr<recob::PFParticle>> pfpVect;
+    // --- Get the candidate slices
+    for(size_t ix=0; ix<sliceVect.size(); ix++){
+        ResetSliceInfo();
+        // --- Get the slice
+        auto & slice = sliceVect[ix];
+
+        // --- Get the slice nu score and check whether we want to correct for it
+
+        // Now I need to get all the hits associated to this flash and get the timing for all of them
+        // Get the slices PFPs
+        double _sliceMaxNuScore = -9999.;
+        pfpVect = slice_pfp_assns.at(slice.key());
+        for(const art::Ptr<recob::PFParticle> &pfp : pfpVect){
+            if(pfp->IsPrimary() &&( std::abs(pfp->PdgCode())==12 || std::abs(pfp->PdgCode())==14 ) ){
+                const std::vector<art::Ptr<larpandoraobj::PFParticleMetadata>> pfpMetaVec = pfp_to_metadata.at(pfp->Self());
+                for (auto const pfpMeta : pfpMetaVec) {
+                    larpandoraobj::PFParticleMetadata::PropertiesMap propertiesMap = pfpMeta->GetPropertiesMap();
+                    if (propertiesMap.count("NuScore")) _fNuScore = propertiesMap.at("NuScore");
+                    if(_fNuScore>_sliceMaxNuScore) _sliceMaxNuScore = _fNuScore;
+                }
+            }
+            std::vector< art::Ptr<recob::Vertex> > vertexVec = pfp_vertex_assns.at(pfp.key());
+            for(const art::Ptr<recob::Vertex> &ver : vertexVec){
+                geo::Point_t xyz_vertex = ver->position();
+                fRecoVx= xyz_vertex.X();
+                fRecoVy= xyz_vertex.Y();
+                fRecoVz= xyz_vertex.Z();
+            }            
+            // If correct light propagation time
+
+
+            // Get the SP associated to the PFP and then get the hits associated to the SP. ---> Hits associated to the PFP
+            //Get the spacepoints associated to the PFParticle
+            std::vector<art::Ptr<recob::SpacePoint>> PFPSpacePointsVect = pfp_sp_assns.at(pfp.key());
+            //Get the SP Hit assns    
+            art::Handle<std::vector<recob::SpacePoint>> eventSpacePoints;
+            std::vector<art::Ptr<recob::SpacePoint>> eventSpacePointsVect;
+            e.getByLabel(fSpacePointLabel, eventSpacePoints);
+            art::fill_ptr_vector(eventSpacePointsVect, eventSpacePoints);
+            art::FindManyP<recob::Hit> SPToHitAssoc (eventSpacePointsVect, e, fSpacePointLabel);
+            for (const art::Ptr<recob::SpacePoint> &SP: PFPSpacePointsVect){
+                std::vector<art::Ptr<recob::Hit>> SPHit = SPToHitAssoc.at(SP.key());
+                if (SPHit.at(0)->WireID().Plane==2){
+                    fSpacePointX.push_back(SP->position().X());
+                    fSpacePointY.push_back(SP->position().Y());
+                    fSpacePointZ.push_back(SP->position().Z());
+                    fSpacePointIntegral.push_back(SPHit.at(0)->Integral());
+
+                    //Fill Bayrcenter Position
+                    if(SP->position().X() < 0){
+                        fChargeWeightX[0] += SP->position().X() * SPHit.at(0)->Integral();
+                        fChargeWeightY[0] += SP->position().Y() * SPHit.at(0)->Integral();
+                        fChargeWeightZ[0] += SP->position().Z() * SPHit.at(0)->Integral();
+                        fChargeTotalWeight[0] += SPHit.at(0)->Integral();
+                    }
+                    else{
+                        fChargeWeightX[1] += SP->position().X() * SPHit.at(0)->Integral();
+                        fChargeWeightY[1] += SP->position().Y() * SPHit.at(0)->Integral();
+                        fChargeWeightZ[1] += SP->position().Z() * SPHit.at(0)->Integral();
+                        fChargeTotalWeight[1] += SPHit.at(0)->Integral();
+                    }
+                }
+            }
+        }
+        //Fill TPC 0 information
+        fChargeBarycenterX[0] = fChargeWeightX[0]/fChargeTotalWeight[0];
+        fChargeBarycenterY[0] = fChargeWeightY[0]/fChargeTotalWeight[0];
+        fChargeBarycenterZ[0] = fChargeWeightZ[0]/fChargeTotalWeight[0];
+        //Fill TPC 1 information
+        fChargeBarycenterX[1] = fChargeWeightX[1]/fChargeTotalWeight[1];
+        fChargeBarycenterY[1] = fChargeWeightY[1]/fChargeTotalWeight[1];
+        fChargeBarycenterZ[1] = fChargeWeightZ[1]/fChargeTotalWeight[1];
+
+        if(_sliceMaxNuScore<fNuScoreThreshold){
+            ResetSliceInfo();
+            continue; // Skip to the next slice if the nu score is below threshold
+        }
+        this->GetPropagationTimeCorrectionPerChannel();
+        
+        // Get the SPECTDC product required to go to the RWM reference frame
+
+        art::Handle<std::vector<sbnd::timing::DAQTimestamp>> tdcHandle;
+        e.getByLabel(fSPECTDCLabel, tdcHandle);
+        if (!tdcHandle.isValid() || tdcHandle->size() == 0){
+            std::cout << "No SPECTDC products found. Skip this event." << std::endl;
+            return;
+        }
+        else{
+            const std::vector<sbnd::timing::DAQTimestamp> tdc_v(*tdcHandle);
+            for (size_t i=0; i<tdc_v.size(); i++){
+                auto tdc = tdc_v[i];
+                const uint32_t  ch = tdc.Channel();
+                const uint64_t  ts = tdc.Timestamp();
+                if(ch == 2) fRWMTime = ts%uint64_t(1e9);
+                if(ch == 4) fEventTriggerTime = ts%uint64_t(1e9);
+            }
+        }
+        // Get all the OpT0 objects associated to the slice
+        std::vector<art::Ptr<recob::OpFlash>> flashFM;
+        if(fFlashMatchingTool == "OpT0Finder" ){
+            const std::vector< art::Ptr<sbn::OpT0Finder> > slcOpT0Finder = slice_opt0finder_assns.at( slice.key() );
+            if(slcOpT0Finder.size() == 0) continue; 
+            size_t OpT0Idx = HighestOpT0ScoreIdx(slcOpT0Finder);
+            // Get the flash OpT0 association
+            art::FindManyP<recob::OpFlash> opflash_opt0finder_assns(opt0Handle, e, fOpT0FinderModuleLabel);
+            flashFM = opflash_opt0finder_assns.at( slcOpT0Finder[OpT0Idx].key() );
+            if(flashFM.size() > 1){
+                throw art::Exception(art::errors::LogicError) << "There are multiple OpFlash objects associated to the same OpT0Finder object. This is not expected.";
+            }
+            _fFMScore = slcOpT0Finder[OpT0Idx]->score;
+            if(_fFMScore < fFMScoreThreshold){
+                ResetSliceInfo();
+                continue;
+            }
+        }
+        else if(fFlashMatchingTool == "BarycenterFM")
+        {
+            const std::vector< art::Ptr<sbn::TPCPMTBarycenterMatch> > slcTPCPMTBarycenter = slice_tpcpmtbarycentermatching_assns.at( slice.key() );
+            if(slcTPCPMTBarycenter.size() == 0) continue;
+            art::FindManyP<recob::OpFlash> opflash_tpcpmtbarycenter_assns(tpcpmtbarycenterfmHandle, e, fTPCPMTBarycenterFMModuleLabel);
+            size_t BFMIdx = HighestBFMScoreIdx(slcTPCPMTBarycenter);
+            flashFM = opflash_tpcpmtbarycenter_assns.at( slcTPCPMTBarycenter[BFMIdx].key() );
+            if(flashFM.size() > 1){
+                throw art::Exception(art::errors::LogicError) << "There are multiple OpFlash objects associated to the same TPCPMTBarycenterFM object. This is not expected.";
+            }
+            _fFMScore = slcTPCPMTBarycenter[BFMIdx]->score;
+            if(_fFMScore < fFMScoreThreshold){
+                ResetSliceInfo();
+                continue;
+            }
+        }
+        else throw art::Exception(art::errors::LogicError) << " Flash matching tool " <<  fFlashMatchingTool << " not supported ." << std::endl; 
+
+
+        // Get the ophits associated to the flash
+        std::vector<art::Ptr<recob::OpHit>> ophitlist;
+        if(flashFM[0]->XCenter()<0)
+        {
+            ophitlist = flashToOpHitAssns_tpc0.at(flashFM[0].key());
+            _mgr = _mgr_tpc0; // Use the TPC 0 flash finder manager
+        }
+        else
+        {
+            ophitlist = flashToOpHitAssns_tpc1.at(flashFM[0].key());
+            _mgr = _mgr_tpc1; // Use the TPC 1 flash finder manager
+        }
+
+        std::vector<recob::OpHit> newOpHitList;
+        std::vector<recob::OpHit> oldOpHitList;
+        for(const auto& ophit : ophitlist) {
+            oldOpHitList.push_back(*ophit);
+        }
+        // Get the list of the corrected OpHits
+        this->CorrectOpHitTime(ophitlist, newOpHitList);
+        // Create the list of ophit lite to be used in the flash finder
+        ::lightana::LiteOpHitArray_t ophits;
+        this->FillLiteOpHit(newOpHitList, ophits);
+        // Create the flash manager
+        auto const flash_v = _mgr.RecoFlash(ophits);
+        double originalFlashTime = flashFM[0]->Time();
+        double newFlashTime = 0.0;
+        for(const auto& lflash :  flash_v) {
+            // Get Flash Barycenter
+            double Ycenter, Zcenter, Ywidth, Zwidth;
+            _flashgeo->GetFlashLocation(lflash.channel_pe, Ycenter, Zcenter, Ywidth, Zwidth);
+            // Get flasht0
+            double flasht0 = lflash.time;
+            // Refine t0 calculation
+            flasht0 = _flasht0calculator->GetFlashT0(lflash.time, GetAssociatedLiteHits(lflash, ophits));
+            recob::OpFlash flash(flasht0, lflash.time_err, flasht0,
+                                ( flasht0) / 1600., lflash.channel_pe,
+                                0, 0, 1, // this are just default values
+                                100., -1., Ycenter, Ywidth, Zcenter, Zwidth);
+            newFlashTime = flasht0;
+            sbn::CorrectedOpFlashTiming correctedOpFlashTiming;
+            correctedOpFlashTiming.OpFlashT0 = originalFlashTime + fEventTriggerTime/1000 - fRWMTime/1000;
+            correctedOpFlashTiming.UpstreamTime_lightonly = originalFlashTime + fEventTriggerTime/1000 - fRWMTime/1000 - (Zcenter/fSpeedOfLight)/1000;
+            correctedOpFlashTiming.UpstreamTime_tpczcorr = originalFlashTime + fEventTriggerTime/1000 - fRWMTime/1000 - (fRecoVz/fSpeedOfLight)/1000;
+            correctedOpFlashTiming.UpstreamTime_propcorr_tpczcorr = newFlashTime + fEventTriggerTime/1000 - fRWMTime/1000 - (fRecoVz/fSpeedOfLight)/1000;
+            correctedOpFlashTiming.FMScore = _fFMScore;
+            correctedOpFlashTiming.SliceNuScore = _sliceMaxNuScore;
+            correctedOpFlashTimes->emplace_back(std::move(correctedOpFlashTiming));
+        }
+
+        art::Ptr<sbn::CorrectedOpFlashTiming> newCorrectedOpFlashTimingPtr = make_correctedopflashtime_ptr(correctedOpFlashTimes->size()-1);
+        newCorrectedOpFlashTimingSliceAssn->addSingle(slice, newCorrectedOpFlashTimingPtr);
+        newCorrectedOpFlashTimingOpFlashAssn->addSingle(flashFM[0], newCorrectedOpFlashTimingPtr);
+        if(fSaveCorrectionTree){
+            this->FillCorrectionTree(newFlashTime, *flashFM[0], oldOpHitList, newOpHitList);
+        }
+    }
+    if(fSaveCorrectionTree) fTree->Fill();
+    ResetEventVars();
+    e.put(std::move(correctedOpFlashTimes));
+    e.put(std::move(newCorrectedOpFlashTimingSliceAssn));
+    e.put(std::move(newCorrectedOpFlashTimingOpFlashAssn));
+    return;
+}
+
+void sbnd::LightPropagationCorrection::beginJob()
+{
+    if(fSaveCorrectionTree)
+    {
+        fTree = tfs->make<TTree>("PMTWaveformFilteranalyzer", "PMT Waveform Filter Analyzer Tree");
+        fTree->Branch("eventID", &fEvent, "eventID/i");
+        fTree->Branch("runID", &fRun, "runID/i");
+        fTree->Branch("subrunID", &fSubrun, "subrunID/i");
+        fTree->Branch("RWMTime", &fRWMTime);
+        fTree->Branch("EventTriggerTime", &fEventTriggerTime);
+        fTree->Branch("NuScore", &fNuScore);
+        fTree->Branch("FMScore", &fFMScore);
+        fTree->Branch("OpFlashTimeOld", &fOpFlashTimeOld);
+        fTree->Branch("OpFlashTimeNew", &fOpFlashTimeNew);
+        fTree->Branch("OpFlashXCenter", &fOpFlashXCenter);
+        fTree->Branch("OpFlashYCenter", &fOpFlashYCenter);
+        fTree->Branch("OpFlashZCenter", &fOpFlashZCenter);
+        fTree->Branch("OpFlashPE", &fOpFlashPE);
+        fTree->Branch("SliceVx", &fSliceVx);
+        fTree->Branch("SliceVy", &fSliceVy);
+        fTree->Branch("SliceVz", &fSliceVz);
+        fTree->Branch("SliceSPX", &fSliceSPX);
+        fTree->Branch("SliceSPY", &fSliceSPY);
+        fTree->Branch("SliceSPZ", &fSliceSPZ);
+        fTree->Branch("OpHitOldTime", &fOpHitOldTime);
+        fTree->Branch("OpHitNewTime", &fOpHitNewTime);
+        fTree->Branch("OpHitPE", &fOpHitPE);
+        fTree->Branch("OpHitOpCh", &fOpHitOpCh);
+    }
+}
+
+void sbnd::LightPropagationCorrection::endJob()
+{
+}
+
+void sbnd::LightPropagationCorrection::ResetEventVars()
+{
+    if(fSaveCorrectionTree)
+    {
+        fEvent = 0;
+        fRun = 0;
+        fSubrun = 0;
+        _fNuScore = 0.0;
+        fRWMTime=-99999.;
+        fEventTriggerTime=-99999.;
+        fNuScore.clear();
+        fFMScore.clear();
+        fOpFlashTimeOld.clear();
+        fOpFlashTimeNew.clear();
+        fOpFlashXCenter.clear();
+        fOpFlashYCenter.clear();
+        fOpFlashZCenter.clear();
+        fOpHitOldTime.clear();
+        fOpHitNewTime.clear();
+        fOpHitPE.clear();
+        fOpHitOpCh.clear();
+        fOpFlashPE.clear();
+        fSliceVx.clear();
+        fSliceVy.clear();
+        fSliceVz.clear();
+        fSliceSPX.clear();
+        fSliceSPY.clear();
+        fSliceSPZ.clear();
+    }
+}
+
+size_t sbnd::LightPropagationCorrection::HighestOpT0ScoreIdx(const std::vector< art::Ptr<sbn::OpT0Finder> > slcFM)
+{
+    // Gets the idx of the OpT0 object with the highest score 
+    double highestOpT0Score = -99999.0; // Initialize to a negative value
+    size_t highestIdx = 0;
+    for(size_t jx=0; jx<slcFM.size(); jx++){
+        if(slcFM[jx]->score > highestOpT0Score){
+            highestOpT0Score = slcFM[jx]->score;
+            highestIdx = jx;
+        }
+    }
+    return highestIdx;
+}
+
+size_t sbnd::LightPropagationCorrection::HighestBFMScoreIdx(const std::vector< art::Ptr<sbn::TPCPMTBarycenterMatch> > slcFM)
+{
+    // Gets the idx of the OpT0 object with the highest score 
+    double highesBFM0Score = -99999.0; // Initialize to a negative value
+    size_t highestIdx = 0;
+    for(size_t jx=0; jx<slcFM.size(); jx++){
+        if(slcFM[jx]->score > highesBFM0Score){
+            highesBFM0Score = slcFM[jx]->score;
+            highestIdx = jx;
+        }
+    }
+    return highestIdx;
+}
+
+void sbnd::LightPropagationCorrection::ResetSliceInfo()
+{
+    _fNuScore=-99999.0;
+    fRecoVx = -99999.0;
+    fRecoVy = -99999.0;
+    fRecoVz = -99999.0;
+    fSpacePointX.clear();
+    fSpacePointY.clear();
+    fSpacePointZ.clear();
+    fSpacePointIntegral.clear();
+    fTimeCorrectionPerChannel.resize(312, 0.0); // Reset the time correction vector for each channel
+    fChargeBarycenterX.assign(2, 0.0);
+    fChargeBarycenterY.assign(2, 0.0);
+    fChargeBarycenterZ.assign(2, 0.0);
+    fChargeWeightX.assign(2, 0.0);
+    fChargeWeightY.assign(2, 0.0);
+    fChargeWeightZ.assign(2, 0.0);
+    fChargeTotalWeight.assign(2, 0.0);
+}
+
+void sbnd::LightPropagationCorrection::GetPropagationTimeCorrectionPerChannel()
+{   
+    // Implementation
+    for(size_t opdet = 0; opdet < fOpDetID.size(); ++opdet) {
+        double _opDetX = fOpDetX[opdet];
+        double _opDetY = fOpDetY[opdet];
+        double _opDetZ = fOpDetZ[opdet];
+        float minPropTime = 999999999.;
+        for(size_t sp=0; sp<fSpacePointX.size(); sp++)
+        {
+            double dx = fSpacePointX[sp] - _opDetX;
+            double dy = fSpacePointY[sp] - _opDetY;
+            double dz = fSpacePointZ[sp] - _opDetZ;
+            double distanceToOpDet = std::sqrt(dx*dx + dy*dy + dz*dz);
+            //double spToCathode = abs(fSpacePointX[sp]); // Distance from space point to cathode in mm
+            //double cathodeToOpDet = std::sqrt(_opDetX*_opDetX + dy*dy + dz*dz); // Distance from cathode to OpDet in mm
+            //float lightPropTimeVIS = spToCathode/fVGroupVUV + cathodeToOpDet/fVGroupVIS; // Speed
+
+            double cathodeToOpDet = std::sqrt(_opDetX*_opDetX + (dy/2)*(dy/2) + (dz/2)*(dz/2)); // Distance from cathode to OpDet in mm
+            double spToCathode = std::sqrt( fSpacePointX[sp]*fSpacePointX[sp] + (dy/2)*(dy/2) + (dz/2)*(dz/2)); // Distance from space point to cathode in mm
+
+            float lightPropTimeVIS = spToCathode/fVGroupVUV + cathodeToOpDet/fVGroupVIS; // Speed
+            float lightPropTimeVUV = distanceToOpDet / fVGroupVUV; // Speed of light in mm/ns for VUV
+            float lightPropTime = std::min(lightPropTimeVIS, lightPropTimeVUV);
+            float partPropTime = std::sqrt((fSpacePointX[sp]-fRecoVx)*(fSpacePointX[sp]-fRecoVx) + (fSpacePointY[sp]-fRecoVy)*(fSpacePointY[sp]-fRecoVy) + (fSpacePointZ[sp]-fRecoVz)*(fSpacePointZ[sp]-fRecoVz))/fSpeedOfLight;
+            float PropTime = lightPropTime + partPropTime;
+            if(PropTime < minPropTime) minPropTime = PropTime;
+        }
+        fTimeCorrectionPerChannel[opdet] = -minPropTime;
+    }
+}
+
+
+void sbnd::LightPropagationCorrection::CorrectOpHitTime(std::vector<art::Ptr<recob::OpHit>> OldOpHitList, std::vector<recob::OpHit> & newOpHitList)
+{
+        int _nophits = 0;
+        _nophits += OldOpHitList.size();
+        for (int i = 0; i < _nophits; ++i) {
+            int opCh = OldOpHitList.at(i)->OpChannel();
+            double channelCorrection = fTimeCorrectionPerChannel[opCh];
+            double newPeakTime = OldOpHitList.at(i)->StartTime()+OldOpHitList.at(i)->RiseTime() + channelCorrection/1000;
+            double newPeakTimeAbs = OldOpHitList.at(i)->PeakTimeAbs()+ channelCorrection/1000;
+            double newStartTime = OldOpHitList.at(i)->StartTime() + channelCorrection/1000;
+            double riseTime = OldOpHitList.at(i)->RiseTime();
+            unsigned int frame = OldOpHitList.at(i)->Frame();
+            double width = OldOpHitList.at(i)->Width();
+            double area = OldOpHitList.at(i)->Area();
+            double amplitude = OldOpHitList.at(i)->Amplitude();
+            double pe = OldOpHitList.at(i)->PE();
+            recob::OpHit newOpHit = recob::OpHit(opCh, newPeakTime, newPeakTimeAbs, newStartTime, riseTime, frame, width, area, amplitude, pe, 0.0);
+            newOpHitList.push_back(newOpHit);
+        }
+}
+
+
+void sbnd::LightPropagationCorrection::FillLiteOpHit(std::vector<recob::OpHit> const& OpHitList, std::vector<::lightana::LiteOpHit_t>& LiteOpHitList)
+{
+    for(auto const& oph : OpHitList) {
+        ::lightana::LiteOpHit_t loph;
+        loph.peak_time = oph.StartTime()+oph.RiseTime();
+        loph.pe = oph.PE();
+        loph.channel = oph.OpChannel();
+        LiteOpHitList.emplace_back(std::move(loph));
+    }
+}
+
+
+void sbnd::LightPropagationCorrection::FillCorrectionTree(double & newFlashTime, recob::OpFlash const& flash, std::vector<recob::OpHit> const& oldOpHitList, std::vector<recob::OpHit> const& newOpHitList){
+
+    fSliceSPX.push_back({});
+    fSliceSPY.push_back({});
+    fSliceSPZ.push_back({});
+    fOpHitOldTime.push_back({});
+    fOpHitNewTime.push_back({});
+    fOpHitPE.push_back({});
+    fOpHitOpCh.push_back({});
+
+    if(fDebug)
+    {
+        for(size_t i=0; i<fSpacePointX.size(); i++){
+            fSliceSPX.back().push_back(fSpacePointX[i]);
+            fSliceSPY.back().push_back(fSpacePointY[i]);
+            fSliceSPZ.back().push_back(fSpacePointZ[i]);
+        }
+        for(size_t i=0; i<oldOpHitList.size(); i++){
+            fOpHitOldTime.back().push_back(oldOpHitList[i].StartTime()+oldOpHitList[i].RiseTime());
+            fOpHitNewTime.back().push_back(newOpHitList[i].StartTime()+newOpHitList[i].RiseTime());
+            fOpHitPE.back().push_back(oldOpHitList[i].PE());
+            fOpHitOpCh.back().push_back(oldOpHitList[i].OpChannel());
+        }
+    }
+
+    fNuScore.push_back(_fNuScore);
+    fFMScore.push_back(_fFMScore);
+    fOpFlashTimeOld.push_back(flash.Time());
+    fOpFlashTimeNew.push_back(newFlashTime);
+    fOpFlashXCenter.push_back(flash.XCenter());
+    fOpFlashYCenter.push_back(flash.YCenter());
+    fOpFlashZCenter.push_back(flash.ZCenter());
+    fOpFlashPE.push_back(std::accumulate(flash.PEs().begin(), flash.PEs().end(), 0.0));
+    fSliceVx.push_back(fRecoVx);
+    fSliceVy.push_back(fRecoVy);
+    fSliceVz.push_back(fRecoVz);
+}
+
+
+::lightana::LiteOpHitArray_t sbnd::LightPropagationCorrection::GetAssociatedLiteHits(::lightana::LiteOpFlash_t lite_flash, ::lightana::LiteOpHitArray_t lite_hits_v)
+  {
+    ::lightana::LiteOpHitArray_t flash_hits_v;
+
+    for(auto const& hitidx : lite_flash.asshit_idx) {
+      flash_hits_v.emplace_back(std::move(lite_hits_v.at(hitidx)));
+    }
+
+    return flash_hits_v;
+  }
\ No newline at end of file
diff --git a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh
new file mode 100644
index 000000000..e5427282a
--- /dev/null
+++ b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh
@@ -0,0 +1,243 @@
+
+// art includes
+#include "art/Framework/Core/EDProducer.h"
+#include "art/Framework/Core/ModuleMacros.h"
+#include "art/Framework/Principal/Event.h"
+#include "art/Framework/Principal/Handle.h"
+#include "art/Framework/Principal/Run.h"
+#include "art/Framework/Principal/SubRun.h"
+#include "art/Persistency/Common/PtrMaker.h"
+#include "art/Utilities/make_tool.h"
+#include "canvas/Utilities/InputTag.h"
+#include "fhiclcpp/ParameterSet.h"
+#include "messagefacility/MessageLogger/MessageLogger.h"
+#include "art_root_io/TFileService.h"
+#include "canvas/Persistency/Common/FindManyP.h"
+#include "canvas/Persistency/Common/FindOneP.h"
+
+// ROOT and C++ includes
+#include <TTree.h>
+#include <string.h>
+
+// Services
+#include "lardataalg/DetectorInfo/DetectorClocksData.h"
+#include "lardata/DetectorInfoServices/DetectorClocksService.h"
+#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
+#include "larcore/Geometry/Geometry.h"
+#include "larcore/Geometry/WireReadout.h"
+#include "lardata/Utilities/AssociationUtil.h"
+
+
+// G4 includes
+#include "lardataobj/Simulation/SimChannel.h"
+#include "lardataobj/Simulation/SimPhotons.h"
+#include "lardataobj/Simulation/SimEnergyDeposit.h"
+
+// Reco includes
+// PDS
+#include "lardataobj/RecoBase/OpFlash.h"
+#include "lardataobj/RecoBase/OpHit.h"
+#include "lardataobj/RawData/OpDetWaveform.h"
+// TPC
+#include "lardataobj/RecoBase/Slice.h"
+#include "lardataobj/RecoBase/SpacePoint.h"
+#include "lardataobj/RecoBase/Hit.h"
+#include "lardataobj/RecoBase/PFParticle.h"
+#include "nusimdata/SimulationBase/MCTruth.h"
+#include "nusimdata/SimulationBase/MCNeutrino.h"
+#include "nusimdata/SimulationBase/MCParticle.h"
+#include "lardataobj/Simulation/SimEnergyDeposit.h"
+#include "lardataobj/Simulation/SimChannel.h"
+#include "lardataobj/RawData/RawDigit.h"
+#include "lardataobj/RecoBase/Track.h"
+#include "lardataobj/RecoBase/Shower.h"
+#include "lardataobj/RecoBase/PFParticle.h"
+#include "lardataobj/RecoBase/SpacePoint.h"
+#include "lardataobj/RecoBase/Cluster.h"
+#include "lardataobj/RecoBase/OpFlash.h"
+#include "lardataobj/RecoBase/OpHit.h"
+#include "lardataobj/RecoBase/Hit.h"
+#include "lardataobj/RawData/RawDigit.h"
+#include "lardataobj/RecoBase/Wire.h"
+#include "lardataobj/RecoBase/Vertex.h"
+#include "lardataobj/RecoBase/Slice.h"
+#include "lardataobj/AnalysisBase/Calorimetry.h"
+#include "lardataobj/AnalysisBase/ParticleID.h"
+#include "lardataobj/AnalysisBase/T0.h"
+#include "lardataobj/RecoBase/PFParticleMetadata.h"
+
+// CRT
+#include "sbnobj/SBND/CRT/CRTTrack.hh"
+
+// Cosmic rejection includes
+#include "sbnobj/Common/Reco/OpT0FinderResult.h"
+#include "sbnobj/Common/Reco/SimpleFlashMatchVars.h"
+#include "sbnobj/Common/Reco/CRUMBSResult.h"
+#include "sbnobj/Common/Reco/TPCPMTBarycenterMatch.h"
+#include "sbnobj/Common/Reco/CorrectedOpFlashTiming.h"
+#include "sbnobj/SBND/Timing/DAQTimestamp.hh"
+#include "lardataobj/AnalysisBase/T0.h"
+
+// Geometry and mapping
+#include "larcore/Geometry/WireReadout.h"
+#include "sbndcode/OpDetSim/sbndPDMapAlg.hh"
+#include <iostream>
+#include <fstream>
+
+// Flash finder utilities
+#include "sbndcode/OpDetReco/OpFlash/FlashFinder/FlashFinderManager.h"
+#include "sbndcode/OpDetReco/OpFlash/FlashFinder/FlashFinderFMWKInterface.h"
+#include "sbndcode/OpDetReco/OpFlash/FlashFinder/PECalib.h"
+#include "sbndcode/OpDetReco/OpFlash/FlashTools/FlashGeoBase.hh"
+#include "sbndcode/OpDetReco/OpFlash/FlashTools/FlashT0Base.hh"
+#include "sbndcode/OpDetReco/OpFlash/FlashTools/DriftEstimatorBase.hh"
+
+#define fXFidCut1 1.5
+#define fXFidCut2 190
+#define fYFidCut 190
+#define fZFidCut1 10
+#define fZFidCut2 490
+
+#define xdet_size 1000
+#define ydet_size 1000
+#define zmindet_size -500
+#define zmaxdet_size 1800
+
+#define fDefaulNeutrinoID 99999
+
+
+namespace sbnd {
+  class LightPropagationCorrection;
+
+}
+
+class sbnd::LightPropagationCorrection : public art::EDProducer {
+public:
+    explicit LightPropagationCorrection(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.
+    LightPropagationCorrection(LightPropagationCorrection const&) = delete;
+    LightPropagationCorrection(LightPropagationCorrection&&) = delete;
+    LightPropagationCorrection& operator=(LightPropagationCorrection const&) = delete;
+    LightPropagationCorrection& operator=(LightPropagationCorrection&&) = delete;
+
+
+private:
+
+    // Required functions.
+    void produce(art::Event & e) override;
+    void endJob() override;
+    void beginJob() override;
+
+    // Selected optional functions.
+    void ResetEventVars();
+    void ResetSliceInfo();
+    size_t HighestOpT0ScoreIdx(const std::vector< art::Ptr<sbn::OpT0Finder> >);
+    size_t HighestBFMScoreIdx(const std::vector< art::Ptr<sbn::TPCPMTBarycenterMatch> > );
+
+    void GetPropagationTimeCorrectionPerChannel();
+    void CorrectOpHitTime(std::vector<art::Ptr<recob::OpHit>> , std::vector<recob::OpHit> & );
+    void FillLiteOpHit(std::vector<recob::OpHit> const& , std::vector<::lightana::LiteOpHit_t>& );
+    void FillCorrectionTree(double & , recob::OpFlash const& , std::vector<recob::OpHit> const& , std::vector<recob::OpHit> const& );
+    ::lightana::LiteOpHitArray_t GetAssociatedLiteHits(::lightana::LiteOpFlash_t , ::lightana::LiteOpHitArray_t );
+
+
+    geo::WireReadoutGeom const& fWireReadout = art::ServiceHandle<geo::WireReadout>()->Get();
+
+    //Flash finder manager
+    ::lightana::FlashFinderManager _mgr;
+    ::lightana::FlashFinderManager _mgr_tpc0;
+    ::lightana::FlashFinderManager _mgr_tpc1;
+
+    // Tool for calculating the OpFlash Y and Z centers
+    std::unique_ptr<lightana::FlashGeoBase> _flashgeo;
+
+    // Tool for calculating the OpFlash t0
+    std::unique_ptr<lightana::FlashT0Base> _flasht0calculator;
+
+    //Vector for PMT position
+    std::vector<double> fOpDetID;
+    std::vector<double> fOpDetX;
+    std::vector<double> fOpDetY;
+    std::vector<double> fOpDetZ;
+
+    std::string fReco2Label;
+    std::string fOpT0FinderModuleLabel;
+    std::string fTPCPMTBarycenterFMModuleLabel;
+    std::string fOpFlashLabel_tpc0;
+    std::string fOpFlashLabel_tpc1;
+    std::string fSpacePointLabel;
+    std::string fOpHitsModuleLabel;
+    std::string fOpFlashNewLabel;
+    std::string fSPECTDCLabel;
+    std::string fFlashMatchingTool;
+    
+    bool fSaveCorrectionTree;
+
+    std::vector<double> fTimeCorrectionPerChannel;
+    double fRecoVx = 0.0;
+    double fRecoVy = 0.0;
+    double fRecoVz = 0.0;
+
+    //Space Point Variables
+    std::vector<double> fSpacePointX;
+    std::vector<double> fSpacePointY;
+    std::vector<double> fSpacePointZ;
+    std::vector<double> fSpacePointIntegral;
+
+    //Charge Barycenter 
+    std::vector<double> fChargeBarycenterX{0.,0.};
+    std::vector<double> fChargeBarycenterY{0.,0.};
+    std::vector<double> fChargeBarycenterZ{0.,0.};
+    std::vector<double> fChargeWeightX{0.,0.};
+    std::vector<double> fChargeWeightY{0.,0.};
+    std::vector<double> fChargeWeightZ{0.,0.};
+    std::vector<double> fChargeTotalWeight{0.,0.};
+
+    double fDriftDistance; // Total Drift Distance
+    double fSpeedOfLight; // Speed of light in mm/ns
+    double fVISLightPropTime;
+    double fKinkDistance;
+    double fVGroupVUV_I;
+    double fVGroupVIS;
+    double fVGroupVUV; // Speed of light in vacuum in mm/ns
+
+    double fNuScoreThreshold;
+    double fFMScoreThreshold;
+
+    bool fDebug;
+
+    art::ServiceHandle<art::TFileService> tfs;
+    TTree *fTree;
+
+    int fEvent;
+    int fRun;
+    int fSubrun;
+    double _fNuScore;
+    double _fFMScore;
+    double fEventTriggerTime=-999999.;
+    double fRWMTime=-999999.;
+    std::vector<double> fNuScore;
+    std::vector<double> fFMScore;
+    std::vector<double> fOpFlashTimeOld;
+    std::vector<double> fOpFlashTimeNew;
+    std::vector<double> fOpFlashXCenter;
+    std::vector<double> fOpFlashYCenter;
+    std::vector<double> fOpFlashZCenter;
+    std::vector<double> fOpFlashPE;
+    std::vector<double> fSliceVx;
+    std::vector<double> fSliceVy;
+    std::vector<double> fSliceVz;
+    std::vector<std::vector<double>> fSliceSPX;
+    std::vector<std::vector<double>> fSliceSPY;
+    std::vector<std::vector<double>> fSliceSPZ;
+    std::vector<std::vector<double>> fOpHitOldTime;
+    std::vector<std::vector<double>> fOpHitNewTime;
+    std::vector<std::vector<double>> fOpHitPE;
+    std::vector<std::vector<int>> fOpHitOpCh;
+};
+
+
+DEFINE_ART_MODULE(sbnd::LightPropagationCorrection)
diff --git a/sbndcode/LightPropagationCorrection/job/CMakeLists.txt b/sbndcode/LightPropagationCorrection/job/CMakeLists.txt
new file mode 100644
index 000000000..eca2a52d6
--- /dev/null
+++ b/sbndcode/LightPropagationCorrection/job/CMakeLists.txt
@@ -0,0 +1 @@
+install_fhicl()
\ No newline at end of file
diff --git a/sbndcode/LightPropagationCorrection/job/run_lightpropagationcorrection.fcl b/sbndcode/LightPropagationCorrection/job/run_lightpropagationcorrection.fcl
new file mode 100644
index 000000000..218248c50
--- /dev/null
+++ b/sbndcode/LightPropagationCorrection/job/run_lightpropagationcorrection.fcl
@@ -0,0 +1,44 @@
+#include "services_sbnd.fcl"
+#include "sbnd_lightpropagationcorrection_config.fcl"
+
+process_name: LIGHTPROPAGATIONCORR  # The process name must NOT contain any underscores
+
+source:
+{
+    module_type: RootInput  # Telling art we want a ROOT input
+    maxEvents:   -1
+}
+
+services:
+{
+    TFileService: { fileName: "lightpropagationcorrection.root" }
+    @table::sbnd_services
+}
+
+outputs:
+{
+    out1:
+    {
+        fileName: "%ifb_lightpropcorr.root"
+        module_type: RootOutput
+        dataTier: "reco"
+        #outputCommands: [ 
+        #                    "keep *_*_*_*"
+        #                    , "drop *_pmtdecoder_*_*"
+        #                ]
+    }
+}
+
+physics:
+{
+    producers:
+    {
+        lightpropcorrection: @local::LightPropagationCorrection
+    }
+
+    prod: [ lightpropcorrection ]
+    stream1: [ out1 ]
+
+    trigger_paths: [ prod ]
+    end_paths: [ stream1 ]
+}
diff --git a/sbndcode/LightPropagationCorrection/job/run_lightpropagationcorrection_ana.fcl b/sbndcode/LightPropagationCorrection/job/run_lightpropagationcorrection_ana.fcl
new file mode 100644
index 000000000..52992e020
--- /dev/null
+++ b/sbndcode/LightPropagationCorrection/job/run_lightpropagationcorrection_ana.fcl
@@ -0,0 +1,33 @@
+#include "services_sbnd.fcl"
+#include "simulationservices_sbnd.fcl"
+#include "messages_sbnd.fcl"
+#include "sam_sbnd.fcl"
+#include "sbnd_lightpropagationcorrectionana_config.fcl"
+
+process_name: CorrectedOpFlashAna
+
+services:
+{
+  TFileService:              { fileName: "correctedopflashana_tree.root" }
+  @table::sbnd_basic_services
+  @table::sbnd_simulation_services
+  ParticleInventoryService: @local::sbnd_particleinventoryservice
+  BackTrackerService:       @local::sbnd_backtrackerservice
+  ParticleInventoryService: @local::standard_particleinventoryservice  
+}
+
+source:
+{
+  module_type: RootInput
+  maxEvents:  -1  # Number of events to create
+}
+
+physics:
+{
+  analyzers:
+  {
+    lightpropcorrectionanatree: @local::lightpropagationcorrection_ana
+  }
+  ana: [lightpropcorrectionanatree]
+  end_paths: [ana]
+}
\ No newline at end of file
diff --git a/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrection_config.fcl b/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrection_config.fcl
new file mode 100644
index 000000000..3277a7bb8
--- /dev/null
+++ b/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrection_config.fcl
@@ -0,0 +1,51 @@
+#include "opticalsimparameterisations_sbnd.fcl"
+#include "sbnd_flashalgo.fcl"
+#include "sbnd_flasht0algo.fcl"
+#include "sbnd_flashgeoalgo.fcl"
+#include "sbnd_flashfinder.fcl"
+
+BEGIN_PROLOG
+
+LightPropagationCorrection:
+{
+    module_type: "LightPropagationCorrection"
+    OpFlashLabel_tpc0: "opflashtpc0"
+    OpFlashLabel_tpc1: "opflashtpc1"
+
+    Reco2Label: "pandora"
+    OpHitsModuleLabel: "ophitpmt"
+    SpacePointLabel: "pandora"
+    OpT0FinderModuleLabel: "opt0finder"
+    TPCPMTBarycenterFMModuleLabel: "tpcpmtbarycentermatching"
+    SPECTDCLabel: "tdcdecoder"
+
+    VGroupVUV: @local::sbnd_vuv_timing_parameterization.vuv_vgroup_mean
+    VGroupVIS: @local::sbnd_vis_timing_parameterization.vis_vmean
+    SpeedOfLight: 29.9
+
+    FlashT0Config: @local::FlashT0SelectedChannels
+    FlashGeoConfig: @local::FlashGeoThreshold
+    AlgoConfig_tpc0: @local::SimpleFlashStandard
+    AlgoConfig_tpc1: @local::SimpleFlashStandard
+    FlashFinderAlgo: @local::SBNDSimpleFlash.FlashFinderAlgo
+
+    SaveSPECTDC: true
+    SaveCorrectionTree: false
+
+    NuScoreThreshold: 0.1
+    FMScoreThreshold: 0.001
+
+    FlashMatchingTool: "BarycenterFM"
+}
+
+LightPropagationCorrection.AlgoConfig_tpc0.TPC: 0
+LightPropagationCorrection.AlgoConfig_tpc1.TPC: 1
+
+LightPropagationCorrectionSCE: @local::LightPropagationCorrection
+LightPropagationCorrectionSCE.Reco2Label: "pandoraSCE"
+LightPropagationCorrectionSCE.SpacePointLabel: "pandoraSCE"
+LightPropagationCorrectionSCE.OpT0FinderModuleLabel: "opt0finderSCE"
+LightPropagationCorrectionSCE.TPCPMTBarycenterFMModuleLabel: "tpcpmtbarycentermatchingSCE"
+
+
+END_PROLOG
\ No newline at end of file
diff --git a/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrectionana_config.fcl b/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrectionana_config.fcl
new file mode 100644
index 000000000..2d0ca7013
--- /dev/null
+++ b/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrectionana_config.fcl
@@ -0,0 +1,10 @@
+
+BEGIN_PROLOG
+
+lightpropagationcorrection_ana:
+{
+  module_type: "LightPropagationCorrectionAna"
+  LightPropCorrectionLabel: "lightpropcorrection"
+}
+
+END_PROLOG
\ No newline at end of file