diff --git a/sbndcode/EvtsFiltering/CMakeLists.txt b/sbndcode/EvtsFiltering/CMakeLists.txt
new file mode 100644
index 000000000..baaf2ccf7
--- /dev/null
+++ b/sbndcode/EvtsFiltering/CMakeLists.txt
@@ -0,0 +1,54 @@
+
+cet_build_plugin(FilterEventID art::module 
+        lardataobj::RecoBase
+        ROOT::Core
+        ROOT::Tree
+        art::Framework_Core
+        art::Framework_Principal
+        art::Utilities
+        art_root_io::tfile_support
+        art_root_io::TFileService_service
+        canvas::canvas
+        fhiclcpp::fhiclcpp
+        cetlib::cetlib
+        cetlib_except::cetlib_except
+        messagefacility::MF_MessageLogger
+        lardataobj::RawData
+        larcore::Geometry_Geometry_service
+        fhiclcpp::fhiclcpp
+        cetlib::cetlib
+        CLHEP::Random
+        ROOT::Geom
+        ROOT::XMLIO
+        ROOT::Gdml
+        ROOT::FFTW
+)
+
+cet_build_plugin(NoFilter art::module 
+        lardataobj::RecoBase
+        ROOT::Core
+        ROOT::Tree
+        art::Framework_Core
+        art::Framework_Principal
+        art::Utilities
+        art_root_io::tfile_support
+        art_root_io::TFileService_service
+        canvas::canvas
+        fhiclcpp::fhiclcpp
+        cetlib::cetlib
+        cetlib_except::cetlib_except
+        messagefacility::MF_MessageLogger
+        lardataobj::RawData
+        larcore::Geometry_Geometry_service
+        fhiclcpp::fhiclcpp
+        cetlib::cetlib
+        CLHEP::Random
+        ROOT::Geom
+        ROOT::XMLIO
+        ROOT::Gdml
+        ROOT::FFTW
+)
+
+install_headers()
+install_source()
+install_fhicl()
diff --git a/sbndcode/EvtsFiltering/FilterEventID_module.cc b/sbndcode/EvtsFiltering/FilterEventID_module.cc
new file mode 100644
index 000000000..a252fbdd0
--- /dev/null
+++ b/sbndcode/EvtsFiltering/FilterEventID_module.cc
@@ -0,0 +1,101 @@
+#include <iostream>
+#include <stdlib.h>
+#include <string>
+#include <vector>
+#include <numeric>
+#include <getopt.h>
+#include <chrono>
+#include <float.h>
+
+
+//some ROOT includes
+#include "TInterpreter.h"
+#include "TROOT.h"
+#include "TH1F.h"
+#include "TH2D.h"
+#include "TCanvas.h"
+#include "TTree.h"
+#include "TFile.h"
+#include "TStyle.h"
+#include "TSystem.h"
+#include "TGraph.h"
+#include "TFFTReal.h"
+#include "TTimeStamp.h"
+#include "TLatex.h"
+#include "TLine.h"
+
+
+// art includes
+#include "canvas/Utilities/InputTag.h"
+#include "canvas/Persistency/Common/FindOneP.h"
+#include "canvas/Persistency/Common/FindManyP.h"
+#include "canvas/Persistency/Common/Ptr.h"
+#include "art/Framework/Core/EDFilter.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 "fhiclcpp/ParameterSet.h" 
+#include "messagefacility/MessageLogger/MessageLogger.h"  
+#include "lardataobj/RecoBase/Hit.h"
+#include "lardataobj/RawData/RawDigit.h"
+#include "lardataobj/RawData/RDTimeStamp.h"
+#include "lardataobj/RawData/raw.h"
+#include "larcore/Geometry/Geometry.h"
+#include "canvas/Persistency/Provenance/Timestamp.h"
+
+#include "art_root_io/TFileService.h"
+
+#include "sbndcode/ChannelMaps/TPC/TPCChannelMapService.h"
+
+
+namespace filt{
+
+  class FilterEventID : public art::EDFilter {
+    public:
+
+      explicit FilterEventID(fhicl::ParameterSet const & pset);
+      virtual bool filter(art::Event& evt) override;
+      void reconfigure(fhicl::ParameterSet const& pset);
+      virtual void beginJob() override;
+
+    private:
+
+      int run;
+      int subrun;
+      int event;
+      std::vector<int> filter_evts;
+    
+  };
+
+
+  FilterEventID::FilterEventID(fhicl::ParameterSet const & pset)
+  : EDFilter(pset)
+  {
+    this->reconfigure(pset);
+    filter_evts = pset.get<std::vector<int> >("filterevts");
+  }
+
+  void FilterEventID::reconfigure(fhicl::ParameterSet const& pset){
+  }
+
+
+  bool FilterEventID::filter(art::Event & evt){
+
+    run = evt.id().run();
+    subrun = evt.id().subRun();
+    event = evt.id().event();
+
+    bool is_in_list = std::find(filter_evts.begin(), filter_evts.end(), event) != filter_evts.end();
+    std::cout << "event: " << event << " is_in_list: " << is_in_list << std::endl;
+    return is_in_list;
+    
+  }
+
+  void FilterEventID::beginJob() {
+  }
+
+  DEFINE_ART_MODULE(FilterEventID)
+
+}
diff --git a/sbndcode/EvtsFiltering/NoFilter_module.cc b/sbndcode/EvtsFiltering/NoFilter_module.cc
new file mode 100644
index 000000000..535044fe2
--- /dev/null
+++ b/sbndcode/EvtsFiltering/NoFilter_module.cc
@@ -0,0 +1,87 @@
+#include <iostream>
+#include <stdlib.h>
+#include <string>
+#include <vector>
+#include <numeric>
+#include <getopt.h>
+#include <chrono>
+#include <float.h>
+
+
+//some ROOT includes
+#include "TInterpreter.h"
+#include "TROOT.h"
+#include "TH1F.h"
+#include "TH2D.h"
+#include "TCanvas.h"
+#include "TTree.h"
+#include "TFile.h"
+#include "TStyle.h"
+#include "TSystem.h"
+#include "TGraph.h"
+#include "TFFTReal.h"
+#include "TTimeStamp.h"
+#include "TLatex.h"
+#include "TLine.h"
+
+
+// art includes
+#include "canvas/Utilities/InputTag.h"
+#include "canvas/Persistency/Common/FindOneP.h"
+#include "canvas/Persistency/Common/FindManyP.h"
+#include "canvas/Persistency/Common/Ptr.h"
+#include "art/Framework/Core/EDFilter.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 "fhiclcpp/ParameterSet.h" 
+#include "messagefacility/MessageLogger/MessageLogger.h"  
+#include "lardataobj/RecoBase/Hit.h"
+#include "lardataobj/RawData/RawDigit.h"
+#include "lardataobj/RawData/RDTimeStamp.h"
+#include "lardataobj/RawData/raw.h"
+#include "larcore/Geometry/Geometry.h"
+#include "canvas/Persistency/Provenance/Timestamp.h"
+
+#include "art_root_io/TFileService.h"
+
+#include "sbndcode/ChannelMaps/TPC/TPCChannelMapService.h"
+
+
+namespace filt{
+
+  class NoFilter : public art::EDFilter {
+    public:
+
+      explicit NoFilter(fhicl::ParameterSet const & pset);
+      virtual bool filter(art::Event& evt) override;
+      void reconfigure(fhicl::ParameterSet const& pset);
+      virtual void beginJob() override;
+
+    private:
+    
+  };
+
+
+  NoFilter::NoFilter(fhicl::ParameterSet const & pset)
+  : EDFilter(pset)
+  {
+    this->reconfigure(pset);
+  }
+
+  void NoFilter::reconfigure(fhicl::ParameterSet const& pset){
+  }
+
+
+  bool NoFilter::filter(art::Event & evt){
+    return true;
+  }
+
+  void NoFilter::beginJob() {
+  }
+
+  DEFINE_ART_MODULE(NoFilter)
+
+}
diff --git a/sbndcode/EvtsFiltering/filtereventid.fcl b/sbndcode/EvtsFiltering/filtereventid.fcl
new file mode 100644
index 000000000..9ee86a2ca
--- /dev/null
+++ b/sbndcode/EvtsFiltering/filtereventid.fcl
@@ -0,0 +1,51 @@
+#include "services_sbnd.fcl"
+
+process_name: FilterEventID
+
+services:
+{
+  TFileService: { fileName: "reco_hist.root"}
+  MemoryTracker:     {}
+  TimeTracker:       {}
+  RandomNumberGenerator: {} 
+  message:              @local::sbnd_message_services_prod_debug
+  FileCatalogMetadata:  @local::art_file_catalog_mc
+  GeometryConfigurationWriter:    {}
+  @table::sbnd_basic_services
+  @table::sbnd_random_services
+}
+
+source:
+{
+  module_type: RootInput
+  maxEvents:  -1
+}
+
+outputs:
+{
+ out1:
+ {
+   module_type: RootOutput
+   fileName:    "%ifb_eventidfiltered.root"
+   dataTier:    "filtered"
+   fastCloning: true
+   SelectEvents: [ filter ]
+ }
+}
+
+physics:
+{
+  filters:
+  {
+    filtereventid: {
+      module_type: FilterEventID
+      filterevts: [2223, 2295]
+    }
+  }
+
+  filter: [ filtereventid ]
+  stream1:  [ out1 ]
+  trigger_paths: [ filter ]
+  end_paths:     [ stream1 ]
+}
+
diff --git a/sbndcode/EvtsFiltering/nofilter.fcl b/sbndcode/EvtsFiltering/nofilter.fcl
new file mode 100644
index 000000000..17bb9e02b
--- /dev/null
+++ b/sbndcode/EvtsFiltering/nofilter.fcl
@@ -0,0 +1,50 @@
+#include "services_sbnd.fcl"
+
+process_name: NoFilter
+
+services:
+{
+  TFileService: { fileName: "reco_hist.root"}
+  MemoryTracker:     {}
+  TimeTracker:       {}
+  RandomNumberGenerator: {} 
+  message:              @local::sbnd_message_services_prod_debug
+  FileCatalogMetadata:  @local::art_file_catalog_mc
+  GeometryConfigurationWriter:    {}
+  @table::sbnd_basic_services
+  @table::sbnd_random_services
+}
+
+source:
+{
+  module_type: RootInput
+  maxEvents:  -1
+}
+
+outputs:
+{
+ out1:
+ {
+   module_type: RootOutput
+   fileName:    "merged.root"
+   dataTier:    "filtered"
+   fastCloning: true
+   SelectEvents: [ filter ]
+ }
+}
+
+physics:
+{
+  filters:
+  {
+    nofilter: {
+      module_type: NoFilter
+    }
+  }
+
+  filter: [ nofilter ]
+  stream1:  [ out1 ]
+  trigger_paths: [ filter ]
+  end_paths:     [ stream1 ]
+}
+
diff --git a/sbndcode/OpDetReco/OpFlash/FlashTools/DriftEstimatorPMTRatio_tool.cc b/sbndcode/OpDetReco/OpFlash/FlashTools/DriftEstimatorPMTRatio_tool.cc
index cf932338d..cf4ff937c 100644
--- a/sbndcode/OpDetReco/OpFlash/FlashTools/DriftEstimatorPMTRatio_tool.cc
+++ b/sbndcode/OpDetReco/OpFlash/FlashTools/DriftEstimatorPMTRatio_tool.cc
@@ -218,7 +218,7 @@ namespace lightana
         drift_distance=Interpolate(pmtratio);
       return drift_distance;
     }
-    else return fDriftCal[fNCalBins-1]; 
+    else return -999999.; 
   }
 
   double DriftEstimatorPMTRatio::GetPropagationTime(double drift){
diff --git a/sbndcode/OpDetReco/OpFlash/FlashTools/FlashGeoThreshold_tool.cc b/sbndcode/OpDetReco/OpFlash/FlashTools/FlashGeoThreshold_tool.cc
index 3d5190370..05fee2eb6 100644
--- a/sbndcode/OpDetReco/OpFlash/FlashTools/FlashGeoThreshold_tool.cc
+++ b/sbndcode/OpDetReco/OpFlash/FlashTools/FlashGeoThreshold_tool.cc
@@ -31,9 +31,15 @@ namespace lightana{
     //Configuration parameters
      struct Config {
 
-       fhicl::Atom<float> Threshold {
-         fhicl::Name("Threshold"),
-         fhicl::Comment("Use channels abnove this thershold")
+       fhicl::Atom<float> ThresholdY {
+         fhicl::Name("ThresholdY"),
+         fhicl::Comment("Use channels above this threshold (Y coord)")
+       };
+
+
+       fhicl::Atom<float> ThresholdZ {
+         fhicl::Name("ThresholdZ"),
+         fhicl::Comment("Use channels above this threshold (Z coord)")
        };
 
        fhicl::Sequence<std::string> PDTypes {
@@ -65,11 +71,13 @@ namespace lightana{
   private:
 
     void ResetVars();
-    void GetCenter(std::map<int, double> PEAcc, double& center, double& width);
+    void GetCenter(std::map<int, double> PEAcc, double& center, double& width, double& threshold);
 
     // Fhicl configuration parameters
     std::vector<std::string> fPDTypes;
-    float fThreshold;
+    double fThresholdY;
+    double fThresholdZ;
+
     bool fNormalizeByPDType;
     unsigned int fWeightExp;
 
@@ -77,7 +85,7 @@ namespace lightana{
     std::map<int, double> fYPEAccumulator;
     std::map<int, double> fZPEAccumulator;
 
-    // Store PD type ar each YZ position
+    // Store PD type at each YZ position
     std::map<int, std::string> fPDTypeByY;
     std::map<int, std::string> fPDTypeByZ;
 
@@ -85,6 +93,12 @@ namespace lightana{
     std::map<std::string, double> fNormFactorsY;
     std::map<std::string, double> fNormFactorsZ;
 
+    // Store the number of ON PMTs at each YZ position
+    std::map<int, unsigned int> fNumONPMTsY_tpc0;
+    std::map<int, unsigned int> fNumONPMTsZ_tpc0;
+    std::map<int, unsigned int> fNumONPMTsY_tpc1;
+    std::map<int, unsigned int> fNumONPMTsZ_tpc1;
+
     // PDS mapping
     opdet::sbndPDMapAlg fPDSMap;
     
@@ -98,29 +112,48 @@ namespace lightana{
 
   FlashGeoThreshold::FlashGeoThreshold(art::ToolConfigTable<Config> const& config)
     : fPDTypes{ config().PDTypes() }
-    , fThreshold{ config().Threshold() }
+    , fThresholdY{ config().ThresholdY() }
+    , fThresholdZ{ config().ThresholdZ() }
     , fNormalizeByPDType{ config().NormalizeByPDType() }
     , fWeightExp{ config().WeightExp() }
   {
 
     // Initialize YZ map
     for(size_t opch=0; opch<::lightana::NOpDets(); opch++){
-      
-      if( std::find(fPDTypes.begin(), fPDTypes.end(), fPDSMap.pdType(opch)) == fPDTypes.end() ) continue;
-      
       ::lightana::OpDetCenterFromOpChannel(opch, fPDxyz);
-      
+      if( std::find(fPDTypes.begin(), fPDTypes.end(), fPDSMap.pdType(opch)) == fPDTypes.end() ) continue;
       fYPEAccumulator[ (int) fPDxyz[1] ] = 0;
       fZPEAccumulator[ (int) fPDxyz[2] ] = 0;
-      
+
+      fNumONPMTsY_tpc0[(int) fPDxyz[1]] = 0;
+      fNumONPMTsY_tpc1[(int) fPDxyz[1]] = 0;
+      fNumONPMTsZ_tpc0[(int) fPDxyz[2]] = 0;
+      fNumONPMTsZ_tpc1[(int) fPDxyz[2]] = 0;
+
       if(fNormalizeByPDType){
+        if( std::find(fPDTypes.begin(), fPDTypes.end(), fPDSMap.pdType(opch)) == fPDTypes.end() ) continue;
         fPDTypeByY[ (int) fPDxyz[1] ] = fPDSMap.pdType(opch);
         fPDTypeByZ[ (int) fPDxyz[2] ] = fPDSMap.pdType(opch);
       }
+    }
 
+    // Initialize accumulators for ON PMTs
+    for(size_t opch=0; opch<::lightana::NOpDets(); opch++){
+      ::lightana::OpDetCenterFromOpChannel(opch, fPDxyz);
+      std::vector<int> SkipChannelList = {39, 66, 67, 71, 85, 86, 87, 92, 115, 138, 141, 170, 197, 217, 218, 221, 222, 223, 226, 245, 248, 249, 302};
+      if(std::find(SkipChannelList.begin(), SkipChannelList.end(), (int) opch) != SkipChannelList.end()) continue; // Skip channels in the list
+      if(fPDxyz[0]<0 && (fPDSMap.pdType(opch)=="pmt_coated" || fPDSMap.pdType(opch)=="pmt_uncoated")) {
+        fNumONPMTsY_tpc0[(int) fPDxyz[1]]++;
+        fNumONPMTsZ_tpc0[(int) fPDxyz[2]]++;
+      }
+      else if(fPDxyz[0]>0 && (fPDSMap.pdType(opch)=="pmt_coated" || fPDSMap.pdType(opch)=="pmt_uncoated")) {
+        fNumONPMTsY_tpc1[(int) fPDxyz[1]]++;
+        fNumONPMTsZ_tpc1[(int) fPDxyz[2]]++;
+      }  
     }
 
-    // Initialize normalixation factors by PD type
+
+    // Initialize normalization factors by PD type
     if(fNormalizeByPDType){
       for(auto & pd:fPDTypes) {
         fNormFactorsY[pd]=0.;
@@ -135,19 +168,24 @@ namespace lightana{
                                             double& Ycenter, double& Zcenter,
                                             double& Ywidth, double& Zwidth)
   {
-
     // Reset variables
     ResetVars();
     Ycenter = Zcenter = fDefCenterValue;
     Ywidth  = Zwidth  = fDefCenterValue;
-    
+
     // Fill PE accumulators
     for (unsigned int opch = 0; opch < pePerOpChannel.size(); opch++) {
       if( std::find(fPDTypes.begin(), fPDTypes.end(), fPDSMap.pdType(opch)) == fPDTypes.end() ) continue;
       // Get physical detector location for this opChannel
       ::lightana::OpDetCenterFromOpChannel(opch, fPDxyz);
-      fYPEAccumulator[(int) fPDxyz[1] ]+=pePerOpChannel[opch];
-      fZPEAccumulator[(int) fPDxyz[2] ]+=pePerOpChannel[opch];
+      if(opch%2==0){
+        fYPEAccumulator[(int) fPDxyz[1] ]+=pePerOpChannel[opch]/fNumONPMTsY_tpc0[(int) fPDxyz[1]];
+        fZPEAccumulator[(int) fPDxyz[2] ]+=pePerOpChannel[opch]/fNumONPMTsZ_tpc0[(int) fPDxyz[2]];
+      }
+      else{
+        fYPEAccumulator[(int) fPDxyz[1] ]+=pePerOpChannel[opch]/fNumONPMTsY_tpc1[(int) fPDxyz[1]];
+        fZPEAccumulator[(int) fPDxyz[2] ]+=pePerOpChannel[opch]/fNumONPMTsZ_tpc1[(int) fPDxyz[2]];
+      }
     }
 
     // Normalize PE accumulators
@@ -161,17 +199,19 @@ namespace lightana{
 
       // Get normalization values for each PD type (Y)
       for(auto & y: fYPEAccumulator){
-        if(y.second>fNormFactorsY[ fPDTypeByY[y.first] ]) 
+        if(y.second>fNormFactorsY[ fPDTypeByY[y.first] ]) {
           fNormFactorsY[ fPDTypeByY[y.first] ] = y.second;
+        }
       }
 
       for(auto & y: fYPEAccumulator)
-        fYPEAccumulator[y.first] = y.second/ fNormFactorsY[ fPDTypeByY[y.first] ]; 
-      
+        fYPEAccumulator[y.first] = y.second/ fNormFactorsY[ fPDTypeByY[y.first] ];
+
       for(auto & z: fZPEAccumulator)
-        fZPEAccumulator[z.first] = z.second/ fNormFactorsZ[ fPDTypeByZ[z.first] ]; 
+        fZPEAccumulator[z.first] = z.second/ fNormFactorsZ[ fPDTypeByZ[z.first] ];
 
     }
+
     else{
 
       double normY = std::max_element(fYPEAccumulator.begin(), fYPEAccumulator.end(), 
@@ -187,9 +227,32 @@ namespace lightana{
  
     }
 
+    double maxVal = 0.0;
+    for (const auto& pair : fYPEAccumulator) {
+        if (pair.second > maxVal)
+            maxVal = pair.second;
+    }
+    if (maxVal > 0.0) {
+      for (auto& pair : fYPEAccumulator) {
+            pair.second /= maxVal;
+      }
+    }
+
+    maxVal = 0.0;
+    for (const auto& pair : fZPEAccumulator) {
+        if (pair.second > maxVal)
+            maxVal = pair.second;
+    }
+    if (maxVal > 0.0) {
+      for (auto& pair : fZPEAccumulator) {
+            pair.second /= maxVal;
+      }
+    }
+
+
     // Get YZ position of selected channels (above threshold)
-    GetCenter(fYPEAccumulator, Ycenter, Ywidth);
-    GetCenter(fZPEAccumulator, Zcenter, Zwidth);
+    GetCenter(fYPEAccumulator, Ycenter, Ywidth, fThresholdY);
+    GetCenter(fZPEAccumulator, Zcenter, Zwidth, fThresholdZ);
 
   }
 
@@ -204,7 +267,7 @@ namespace lightana{
     }
   }
 
-  void FlashGeoThreshold::GetCenter(std::map<int, double> PEAcc, double& center, double& width){
+  void FlashGeoThreshold::GetCenter(std::map<int, double> PEAcc, double& center, double& width, double& threshold){
 
     // set variables
     center = 0.;
@@ -227,7 +290,7 @@ namespace lightana{
       sum2+=weight*pe.first*pe.first;
 
       // For OpFlash center consider only channels above ceertain threshold
-      if(pe.second>fThreshold){  
+      if(pe.second>threshold){  
         center+=weight*pe.first;
         weightNormCenter+=weight;
       }
diff --git a/sbndcode/OpDetReco/OpFlash/job/sbnd_flashgeoalgo.fcl b/sbndcode/OpDetReco/OpFlash/job/sbnd_flashgeoalgo.fcl
index f7385048f..18051a4c2 100644
--- a/sbndcode/OpDetReco/OpFlash/job/sbnd_flashgeoalgo.fcl
+++ b/sbndcode/OpDetReco/OpFlash/job/sbnd_flashgeoalgo.fcl
@@ -10,7 +10,8 @@ FlashGeoThreshold:
 {
   tool_type: "FlashGeoThreshold"
   PDTypes: ["pmt_coated", "pmt_uncoated"]
-  Threshold: 0.8
+  ThresholdY: 0.75
+  ThresholdZ: 0.45
   NormalizeByPDType: true
   WeightExp: 2
 }
diff --git a/sbndcode/TPCPMTBarycenterMatching/CMakeLists.txt b/sbndcode/TPCPMTBarycenterMatching/CMakeLists.txt
index 7158397dd..9a409d8f0 100644
--- a/sbndcode/TPCPMTBarycenterMatching/CMakeLists.txt
+++ b/sbndcode/TPCPMTBarycenterMatching/CMakeLists.txt
@@ -1,5 +1,6 @@
 set(	MODULE_LIBRARIES
  			sbnobj::Common_Trigger
+			sbndcode_OpDetSim
             larcorealg::Geometry
 			larcore::Geometry_Geometry_service
 			lardataobj::RecoBase
@@ -7,6 +8,8 @@ set(	MODULE_LIBRARIES
 			larcoreobj::SummaryData
 			larana::OpticalDetector_OpHitFinder
 			larreco::Calorimetry
+			larsim::PhotonPropagation
+            larsim::PhotonPropagation_PhotonVisibilityService_service
 			lardata::Utilities
 			nusimdata::SimulationBase
 			nurandom::RandomUtils_NuRandomService_service
@@ -31,6 +34,7 @@ set(	MODULE_LIBRARIES
 			ROOT::FFTW
 			ROOT::Core
             ROOT::Tree
+			Eigen3::Eigen
 )
 
 cet_build_plugin(TPCPMTBarycenterMatchProducer art::module SOURCE TPCPMTBarycenterMatching_module.cc LIBRARIES ${MODULE_LIBRARIES})
diff --git a/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc b/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc
index 8da5761d5..a9f260fe3 100644
--- a/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc
+++ b/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc
@@ -39,7 +39,10 @@
 #include "larcore/Geometry/Geometry.h"
 #include "larcore/Geometry/WireReadout.h"
 #include "larcorealg/Geometry/GeometryCore.h"
+#include "larsim/PhotonPropagation/SemiAnalyticalModel.h"
 #include "larcore/CoreUtils/ServiceUtil.h"
+#include "sbndcode/OpDetSim/sbndPDMapAlg.hh"
+
 
 //Data product includes
 #include "lardataobj/RecoBase/OpHit.h"
@@ -48,13 +51,19 @@
 #include "lardataobj/RecoBase/Hit.h"
 #include "lardataobj/RecoBase/SpacePoint.h"
 #include "lardataobj/RecoBase/PFParticle.h"
+#include "lardataobj/RecoBase/Cluster.h"
 #include "lardataobj/AnalysisBase/T0.h"
 #include "lardataobj/RawData/TriggerData.h"
+#include "lardataobj/RecoBase/PFParticleMetadata.h"
 #include "sbnobj/Common/Reco/TPCPMTBarycenterMatch.h"
 
 //ROOT includes
+#include <Eigen/Dense>
+#include <vector>
+#include <numeric>
+#include <cmath>
 #include "TTree.h"
-#include "TVector3.h"
+
 
 #include <cmath> // std::hypot(), std::abs(), std::sqrt()
 #include <iostream>
@@ -207,16 +216,31 @@ class TPCPMTBarycenterMatchProducer : public art::EDProducer {
   void updateChargeVars(double sumCharge, TVector3 const& sumPos, TVector3 const& sumPosSqr, std::array<double, 3> const& triggerFlashCenter); ///< Update slice-level data members with charge and trigger match info
   void updateFlashVars(art::Ptr<recob::OpFlash> flash);                      ///< Update slice-level data members with best match info
   void updateMatchInfo(sbn::TPCPMTBarycenterMatch& matchInfo);                                      ///< Update match product with slice-level data members
- 
+  void GetPCA(std::vector<double> const& x, std::vector<double> const& y, std::vector<double> const& weight, std::vector<double>& PCA ); ///< Get the PCA
+  double GetSliceCharge(const std::vector<art::Ptr<recob::Hit>> &tpcHitsVec, const detinfo::DetectorPropertiesData det_prop, int tpc);
+  double GetFlashLight(double flash_pe, std::vector<double>& total_dir_visibility, std::vector<double>& total_ref_visibility);
+  void CreateOpHitList( std::vector<art::Ptr<recob::OpHit>> ophitlist, std::vector<double>& ophit_z, std::vector<double>& ophit_y, std::vector<double>& ophit_weight);
+
   // Input parameters
   std::vector<std::string>  fInputTags;            ///< Suffix added onto fOpFlashLabel and fPandoraLabel, used by ICARUS for separate cryostat labels but could be empty
   std::vector<std::string>  fOpFlashesModuleLabel;         ///< Label for PMT reconstruction products
   std::string               fPandoraLabel;         ///< Label for Pandora output products
+  std::string               fOpT0Label;
   bool                      fCollectionOnly;       ///< Only use TPC spacepoints from the collection plane
+  double                    fDistanceCandidateFlashes; ///< Maximum distance between candidate flashes to be considered for matching (cm)
+  std::vector<double>       fCalAreaConst;         /// Calibration area constants for wire plane
+  std::vector<int>          fSkipChannelList;
+  double                    fOpDetVUVEff;           // Efficiencies for PMT detection
+  double                    fOpDetVISEff;           // Efficiencies for PMT detection
   bool                      fVerbose;              ///< Print extra info
   bool                      fFillMatchTree;        ///< Fill an output TTree in the supplemental file
   bool                      fDo3DMatching;         ///< Wether to perform the matching in 3D or 2D
-
+  std::vector<double>       fLightChargeRatioBounds; ///< Vector to store the distance between the barycenter of the charge and the barycenter of the light for each slice
+  double                    fXError;
+  double                    fYError;
+  double                    fZError;
+  double                    fAngleError;
+  std::vector<double>       fFlashVetoWindow;
   // Event-level data members
   int                       fRun;                  ///< Number of the run being processed
   int                       fEvent;                ///< Number of the event being processed
@@ -232,6 +256,7 @@ class TPCPMTBarycenterMatchProducer : public art::EDProducer {
   double                    fChargeWidthY;         ///< Weighted standard deviation of Y position of spacepoints (cm)
   double                    fChargeWidthZ;         ///< Weighted standard deviation of Z position of spacepoints (cm)
   double                    fFlashTime;            ///< Matched OpFlash time (us)
+  double                    fFlashTimeDebug;         ///< Matched OpFlash time error (us)
   double                    fFlashPEs;             ///< Brightness of matched flash (photoelectrons)
   double                    fFlashCenterX;         ///< Flash position X obtained through eta_pmt curve (cm)
   double                    fFlashCenterY;         ///< Weighted mean Y postion of hit PMTs (cm)
@@ -239,16 +264,31 @@ class TPCPMTBarycenterMatchProducer : public art::EDProducer {
   double                    fFlashWidthY;          ///< Weighted standard deviation of Y postion of hit PMTs (cm)
   double                    fFlashWidthZ;          ///< Weighted standard deviation of Z postion of hit PMTs (cm)
   double                    fDeltaT;               ///< | Matched flash time - charge T0 | when available (us)
-  double                    fDeltaX;               ///< | Matched flash X center - charge X center | (cm)
   double                    fDeltaY;               ///< | Matched flash Y center - charge Y center | (cm)
   double                    fDeltaZ;               ///< | Matched flash Z center - charge Z center | (cm)
   double                    fRadius;               ///< Hypotenuse of DeltaY and DeltaZ *parameter minimized by matching* (cm)
-  double                    fDeltaX_Trigger;       ///< | Triggering flash X center - charge X center | (cm)
+  double                    fChi2;                 ///< Chi2 to be minimized when matching flash to slice (dimensionless)
+  double                    fScore;                ///< Score to be maximized when matching flash to slice (dimensionless)
+  double                    fAngle;                ///< Angle between charge PCA and light PCA (degrees)
   double                    fDeltaY_Trigger;       ///< | Triggering flash Y center - charge Y center | (cm)
   double                    fDeltaZ_Trigger;       ///< | Triggering flash Z center - charge Z center | (cm)
   double                    fRadius_Trigger;       ///< Hypotenuse of DeltaY_Trigger and DeltaZ_Trigger (cm)
-  
   TTree*                    fMatchTree;            ///< Tree to store all match information
+  // Geometry service
+  geo::WireReadoutGeom const& fWireReadout = art::ServiceHandle<geo::WireReadout>()->Get();
+  opdet::sbndPDMapAlg fPDSMap;
+
+  //Vector for PMT position
+  std::vector<double> fOpDetID;
+  std::vector<int> fOpDetType;
+  std::vector<double> fOpDetX;
+  std::vector<double> fOpDetY;
+  std::vector<double> fOpDetZ;
+
+  // Semi-analytical model for VUV and VIS light
+  std::unique_ptr<phot::SemiAnalyticalModel> _semi_model;
+  fhicl::ParameterSet _vuv_params;
+  fhicl::ParameterSet _vis_params;
 
 };
 
@@ -259,9 +299,19 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet
   fOpFlashesModuleLabel(p.get<std::vector<std::string>>("OpFlashesModuleLabel")),
   fPandoraLabel(p.get<std::string>("PandoraLabel")),
   fCollectionOnly(p.get<bool>("CollectionOnly", true)),
-  fVerbose(p.get<bool>("Verbose", false)),
-  fFillMatchTree(p.get<bool>("FillMatchTree", false)),
-  fDo3DMatching(p.get<bool>("Do3DMatching",true))
+  fDistanceCandidateFlashes(p.get<double>("DistanceCandidateFlashes")), // cm
+  fCalAreaConst(p.get<std::vector<double>>("CalAreaConst")),
+  fOpDetVUVEff (p.get<double>("OpDetVUVEff")),
+  fOpDetVISEff (p.get<double>("OpDetVISEff")),
+  fVerbose(p.get<bool>("Verbose")),
+  fFillMatchTree(p.get<bool>("FillMatchTree")),
+  fDo3DMatching(p.get<bool>("Do3DMatching")),
+  fLightChargeRatioBounds(p.get<std::vector<double>>("LightChargeRatioBounds")),
+  fXError(p.get<double>("XError")), // cm
+  fYError(p.get<double>("YError")), // cm
+  fZError(p.get<double>("ZError")), // cm
+  fAngleError(p.get<double>("AngleError")), // deg
+  fFlashVetoWindow(p.get<std::vector<double>>("FlashVetoWindow")) // us
   {
   // Call appropriate produces<>() functions here.
 
@@ -274,14 +324,15 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet
     << "Module has been missconfigured. Number of OpFlashesModuleLabel must be 2!";
   }
 
+  art::ServiceHandle<art::TFileService> tfs;
   if ( fFillMatchTree ) {
-    art::ServiceHandle<art::TFileService> tfs;
-    fMatchTree = tfs->make<TTree>("matchTree","TPC Slice - OpFlash Matching Analysis");
+      art::ServiceHandle<art::TFileService> tfs;
+      fMatchTree = tfs->make<TTree>("matchTree","TPC Slice - OpFlash Matching Analysis");
 
     //Event Info
     fMatchTree->Branch("run",                 &fRun,                 "run/I"                );
     fMatchTree->Branch("event",               &fEvent,               "event/I"              );
-    fMatchTree->Branch("cryo",                &fTPC,                "cryo/I"               );
+    fMatchTree->Branch("cryo",                &fTPC,                "cryo/I"                );
     fMatchTree->Branch("sliceNum",            &fSliceNum,            "sliceNum/I"           );
 
     //Charge Info
@@ -305,17 +356,36 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet
 
     //Match Quality Info
     fMatchTree->Branch("deltaT",              &fDeltaT,              "deltaT/d"             );
-    fMatchTree->Branch("deltaX",              &fDeltaX,              "deltaX/d"             );
     fMatchTree->Branch("deltaY",              &fDeltaY,              "deltaY/d"             );
     fMatchTree->Branch("deltaZ",              &fDeltaZ,              "deltaZ/d"             );
     fMatchTree->Branch("radius",              &fRadius,              "radius/d"             );
-    fMatchTree->Branch("deltaX_Trigger",      &fDeltaX_Trigger,      "deltaX_Trigger/d"     );
+    fMatchTree->Branch("angle",               &fAngle,               "angle/d"              );
+    fMatchTree->Branch("chi2",                &fChi2,                "chi2/d"               );
+    fMatchTree->Branch("score",               &fScore,               "score/d"              );
     fMatchTree->Branch("deltaZ_Trigger",      &fDeltaZ_Trigger,      "deltaZ_Trigger/d"     );
     fMatchTree->Branch("deltaY_Trigger",      &fDeltaY_Trigger,      "deltaY_Trigger/d"     );
-    fMatchTree->Branch("radius_Trigger",      &fRadius_Trigger,      "dadius_Trigger/d"     );
+    fMatchTree->Branch("radius_Trigger",      &fRadius_Trigger,      "radius_Trigger/d"     );
 
   } //End MatchTree
 
+  //Fill the OpDet positions
+  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());
+    if(fPDSMap.pdType(opch)=="pmt_coated") fOpDetType.push_back(0);
+    else if(fPDSMap.pdType(opch)=="pmt_uncoated") fOpDetType.push_back(1);
+    else if(fPDSMap.pdType(opch)=="xarapuca_vuv") fOpDetType.push_back(2);
+    else if(fPDSMap.pdType(opch)=="xarapuca_vis") fOpDetType.push_back(3);
+    else fOpDetType.push_back(-1);
+  }
+
+  _vuv_params = p.get<fhicl::ParameterSet>("VUVHits");
+  _vis_params = p.get<fhicl::ParameterSet>("VIVHits");
+  _semi_model = std::make_unique<phot::SemiAnalyticalModel>(_vuv_params, _vis_params, true, false);
+
 }
 
 void TPCPMTBarycenterMatchProducer::produce(art::Event& e)
@@ -324,93 +394,148 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e)
   fEvent  = e.id().event();
   fRun    = e.run();
 
+  // Detector properties and clocks
+  auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(e);
+  auto const det_prop = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(e, clock_data);
+
+
   //Initialize new data products
   auto matchInfoVector = std::make_unique< std::vector<sbn::TPCPMTBarycenterMatch> >();
   art::PtrMaker< sbn::TPCPMTBarycenterMatch > const makeInfoPtr(e); 
   auto sliceAssns = std::make_unique< art::Assns<sbn::TPCPMTBarycenterMatch, recob::Slice> >();
   auto flashAssns = std::make_unique< art::Assns<sbn::TPCPMTBarycenterMatch, recob::OpFlash> >();
 
-    /* ~~~~~~~~~~~~~~~~~~~~ Flash Section
- *
- * Here we gather the OpFlashes found in this cryostat and their OpHits
- * We iterate through the flashes to identify a triggering flash
- */
-  std::array<double, 3> triggerFlashCenter ;
-  for ( size_t tpc=0; tpc<2 ; tpc ++) {
-    //Fetch the flashes and their associated hits, pointer vector needed for generating associations
-    art::Handle flashHandle = e.getHandle<std::vector<recob::OpFlash>>(fOpFlashesModuleLabel[tpc]);
-    int nFlashes = (*flashHandle).size();
-    triggerFlashCenter = {-9999., -9999., -9999.};
-    double minTime = 99999.;
-    for (int i = 0; i < nFlashes; ++i) {
-        const recob::OpFlash &flash = (*flashHandle)[i];
-        double _flashtime = flash.AbsTime();
-        if ( std::abs(_flashtime)< minTime) {
-            minTime = std::abs(_flashtime);
-            triggerFlashCenter = {flash.XCenter(), flash.YCenter(), flash.ZCenter()};
-        }
+  /* ~~~~~~~~~~~~~~~~~~~~ Flash Section
+  *
+  * Here we gather the OpFlashes found in this cryostat and their OpHits
+  * We iterate through the flashes to identify a triggering flash
+  */
+    std::array<double, 3> triggerFlashCenter;
+    for ( size_t tpc=0; tpc<2 ; tpc ++) {
+      //Fetch the flashes and their associated hits, pointer vector needed for generating associations
+      art::Handle flashHandle = e.getHandle<std::vector<recob::OpFlash>>(fOpFlashesModuleLabel[tpc]);
+      int nFlashes = (*flashHandle).size();
+      triggerFlashCenter = {-9999., -9999., -9999.};
+      double minTime = 99999.;
+      for (int i = 0; i < nFlashes; ++i) {
+          const recob::OpFlash &flash = (*flashHandle)[i];
+          double _flashtime = flash.AbsTime();
+          if ( std::abs(_flashtime)< minTime) {
+              minTime = std::abs(_flashtime);
+              triggerFlashCenter = {flash.XCenter(), flash.YCenter(), flash.ZCenter()};
+          }
+      }
     }
-  }
+
 
   // TODO (acastill): evaluate if there are two flashes in time coincidence and treat them as one.
 
-/* ~~~~~~~~~~~~~~~~~~~~ TPC Section
- * Here we start by gathering the Slices in the event
- * For each slice, the charge centroid is first calculated
- * Then we iterate through flashes to identify the best match flash
- * If a triggering flash was found in this cyrostat, the barycenter distance to the triggering flash is also stored
- */
+  /* ~~~~~~~~~~~~~~~~~~~~ TPC Section
+  * Here we start by gathering the Slices in the event
+  * For each slice, the charge centroid is first calculated
+  * Then we iterate through flashes to identify the best match flash
+  * If a triggering flash was found in this cyrostat, the barycenter distance to the triggering flash is also stored
+  */
+
+  //Fetch slices, TPC hits, and PFPs; pointer vector needed for generating associations
+  ::art::Handle<std::vector<recob::Slice>> sliceHandle;
+  e.getByLabel(fPandoraLabel, sliceHandle);
 
-    //Fetch slices, TPC hits, and PFPs; pointer vector needed for generating associations
-  art::Handle const sliceHandle
-    = e.getHandle<std::vector<recob::Slice>>(fPandoraLabel);
-  art::FindManyP<recob::Hit> fmTPCHits(sliceHandle, e, fPandoraLabel);
-  art::FindManyP<recob::PFParticle> fmPFPs(sliceHandle, e, fPandoraLabel);
+  // Slice to PFP assns
   unsigned nSlices = (*sliceHandle).size();
+  ::art::Handle<std::vector<recob::PFParticle>> pfpHandle;
+  e.getByLabel(fPandoraLabel, pfpHandle);
+
+  // PFP to metadata assns
+  art::FindManyP<recob::PFParticle> slice_pfp_assns(sliceHandle, e, fPandoraLabel);
+
+  // Hit to PFP assns
+  art::FindManyP<recob::Cluster> pfp_cluster_assns(pfpHandle, e, fPandoraLabel);
+  ::art::Handle<std::vector<recob::Cluster>> clusterHandle;
+  e.getByLabel(fPandoraLabel, clusterHandle);
+
+  // Cluster to Hit assns
+  art::FindManyP<recob::Hit> cluster_hit_assns (clusterHandle, e, fPandoraLabel);
+
+    
 
   //For slice...
   for ( unsigned j = 0; j < nSlices; j++ ) {
-  //For each TPC
-  std::vector<sbn::TPCPMTBarycenterMatch> sliceMatchInfoVector;
-  std::vector<art::Ptr<sbn::TPCPMTBarycenterMatch>>  infoPtrVector;
-  std::vector<art::Ptr<recob::OpFlash>> flashPtrVector;
-  fSliceNum = j;
-  const art::Ptr<recob::Slice> slicePtr { sliceHandle, j };
+    //For each TPC
+    std::vector<sbn::TPCPMTBarycenterMatch> sliceMatchInfoVector;
+    std::vector<art::Ptr<sbn::TPCPMTBarycenterMatch>>  infoPtrVector;
+    std::vector<art::Ptr<recob::OpFlash>> flashPtrVector;
+    fSliceNum = j;
+    const art::Ptr<recob::Slice> slicePtr { sliceHandle, j };
+    //Vector for recob PFParticles
+    std::vector<art::Ptr<recob::Hit>> tpcHitsVec;
+    std::vector<art::Ptr<recob::PFParticle>> pfpVect = slice_pfp_assns.at(j);
+
+    // Get the hits associated to the PFParticles in this slice
+    for(const art::Ptr<recob::PFParticle> &pfp : pfpVect){
+      std::vector<art::Ptr<recob::Cluster>> cluster_v = pfp_cluster_assns.at(pfp.key());
+      for(size_t i=0; i<cluster_v.size(); i++){
+        std::vector<art::Ptr<recob::Hit>> hitVect = cluster_hit_assns.at(cluster_v[i].key());
+        tpcHitsVec.insert(tpcHitsVec.end(), hitVect.begin(), hitVect.end());
+      }
+    }
+
+    int nPFPs = pfpVect.size();
+    //Retrieve Pandora's T0 for this slice if available, same for every PFP in slice so we only need one
+    if ( nPFPs != 0 ) {
+      art::FindOne<anab::T0> f1T0( {pfpVect.at(0)}, e, fPandoraLabel);
+      if ( f1T0.at(0).isValid() ) {
+        fChargeT0 = f1T0.at(0).ref().Time() / 1e3;
+      }
+    }
+
     for ( size_t tpc=0; tpc<2 ; tpc ++) {
       fTPC = tpc;
       InitializeSlice();
       sbn::TPCPMTBarycenterMatch sliceMatchInfo;
       updateMatchInfo(sliceMatchInfo);
-
-      const std::vector<art::Ptr<recob::Hit>> &tpcHitsVec = fmTPCHits.at(j);
-      const std::vector<art::Ptr<recob::PFParticle>> &pfpsVec = fmPFPs.at(j);
       art::FindOne<recob::SpacePoint> f1SpacePoint(tpcHitsVec, e, fPandoraLabel);
 
-      int nHits = tpcHitsVec.size();
-      int nPFPs = pfpsVec.size();
+      std::vector<double> hit_z;
+      std::vector<double> hit_y;
+      std::vector<double> hit_weight;
 
-      //Retrieve Pandora's T0 for this slice if available, same for every PFP in slice so we only need one
-      if ( nPFPs != 0 ) {
-        art::FindOne<anab::T0> f1T0( {pfpsVec.at(0)}, e, fPandoraLabel);
-        if ( f1T0.at(0).isValid() ) {
-          fChargeT0 = f1T0.at(0).ref().Time() / 1e3;
-        }
-      }
+      int nHits = tpcHitsVec.size();
 
       double thisCharge;
       double sumCharge = 0.;
       TVector3 sumPos {0.,0.,0.};
       TVector3 sumPosSqr {0.,0.,0.};
 
+      size_t maxChargePlaneIdx=99999;
+
+      // If we are not using collection plane only, we need to find the plane with the most charge
+      if(!fCollectionOnly){
+        std::vector<double> PlaneCharge(3, 0.); // Vector to store the charge for each plane
+        //Loop to get the charge of each plane
+        for ( int k = 0; k < nHits; k++ ) {
+          if ( !f1SpacePoint.at(k).isValid() ) continue;
+          // If hit does not belong to correct tpc then skip
+          const art::Ptr<recob::Hit> &tpcHit = tpcHitsVec.at(k);
+          const recob::SpacePoint point = f1SpacePoint.at(k).ref();
+          TVector3 const thisPoint = point.XYZ();
+          if ((tpc == 0) == (thisPoint.X() > 0)) continue; // Skip if the point is not in the TPC we are considering
+          int plane = tpcHit->WireID().Plane;
+          PlaneCharge[plane] += tpcHit->Integral();
+        }
+        // Idx to select the plane with the most charge
+        maxChargePlaneIdx = std::distance(PlaneCharge.begin(), std::max_element(PlaneCharge.begin(), PlaneCharge.end()));
+      }
+
       //For hit...
       for ( int k = 0; k < nHits; k++ ) {
         // If hit does not belong to correct tpc then skip 
         const art::Ptr<recob::Hit> &tpcHit = tpcHitsVec.at(k);
 
         //Only use hits with associated SpacePoints, and optionally only collection plane hits
-        if ( fCollectionOnly && tpcHit->SignalType() != geo::kCollection ) continue;
         if ( !f1SpacePoint.at(k).isValid() ) continue;
-
+        if ( fCollectionOnly && tpcHit->WireID().Plane != 2 ) continue;
+        else if ( !fCollectionOnly && tpcHit->WireID().Plane != maxChargePlaneIdx) continue; // If not using collection plane only, skip hits not in the plane with the most charge
         const recob::SpacePoint point = f1SpacePoint.at(k).ref();
         TVector3 const thisPoint = point.XYZ();
         if ((tpc == 0) == (thisPoint.X() > 0)) continue; // Skip if the point is not in the TPC we are considering
@@ -419,6 +544,10 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e)
         sumCharge += thisCharge;
         sumPos += thisPoint * thisCharge;
         sumPosSqr += thisPointSqr * thisCharge;
+        // Store hit information for PCA
+        hit_y.push_back(thisPoint.Y());
+        hit_z.push_back(thisPoint.Z());
+        hit_weight.push_back(thisCharge);
       } //End for hit
 
       //No charge found in slice...
@@ -431,27 +560,61 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e)
       //Update charge variables
       updateChargeVars(sumCharge, sumPos, sumPosSqr, triggerFlashCenter);
       updateMatchInfo(sliceMatchInfo);
+      
+      std::vector<double> ChargePCA = {0.,0.};
+      GetPCA(hit_z, hit_y, hit_weight, ChargePCA);
+
+      //Get the charge of the slice 
+      double sliceCharge = this->GetSliceCharge(tpcHitsVec, det_prop, tpc);
+      //Get the visibility map of the charge barycenter
+      geo::Point_t SliceXYZ = {fChargeCenterX, fChargeCenterY, fChargeCenterZ};
+      std::vector<double> direct_visibility;
+      std::vector<double> reflect_visibility;
+      _semi_model->detectedDirectVisibilities(direct_visibility, SliceXYZ);
+      _semi_model->detectedReflectedVisibilities(reflect_visibility, SliceXYZ);
 
       int matchIndex = -1;
       double minDistance = 1e6;
       double thisFlashCenterX, thisFlashCenterY, thisFlashCenterZ, thisDistance;
 
-      art::Handle flashHandle = e.getHandle<std::vector<recob::OpFlash>>(fOpFlashesModuleLabel[tpc]);
-      int nFlashes = (*flashHandle).size();
+      // --- Read Recob OpFlash
+      ::art::Handle<std::vector<recob::OpFlash>> flashHandle;
+      e.getByLabel(fOpFlashesModuleLabel[tpc], flashHandle);
+      std::vector< art::Ptr<recob::OpFlash> > flashVect;
+      art::fill_ptr_vector(flashVect, flashHandle);
+
+      //OpHit OpFlash assns
+      art::FindManyP<recob::OpHit> flash_ophit_assns(flashHandle, e, fOpFlashesModuleLabel[tpc]);
+
+      //Vector to store the idxs of the candidate flashes
+      std::vector<int> candidateFlashIdxs;
+      int nFlashes = flashVect.size();
+
       //For flash...
       for ( int m = 0; m < nFlashes; m++ ) {
-        const recob::OpFlash &flash = (*flashHandle).at(m);
+        auto & flash = flashVect[m];
+
+        double flashTime = flash->Time();
+
+        if(flashTime<fFlashVetoWindow[0] || flashTime>fFlashVetoWindow[1]) continue;
+
         //Find index of flash that minimizes barycenter distance in XYZ place
-        thisFlashCenterX = flash.XCenter();
-        thisFlashCenterY = flash.YCenter();
-        thisFlashCenterZ = flash.ZCenter();
+        thisFlashCenterX = flash->XCenter();
+        thisFlashCenterY = flash->YCenter();
+        thisFlashCenterZ = flash->ZCenter();
+
         if(fDo3DMatching)
           thisDistance = std::hypot( (thisFlashCenterX - fChargeCenterX), (thisFlashCenterY - fChargeCenterY), (thisFlashCenterZ - fChargeCenterZ) );
         else thisDistance = std::hypot( (thisFlashCenterY - fChargeCenterY), (thisFlashCenterZ - fChargeCenterZ) );
+
         if ( thisDistance < minDistance ) {
           minDistance = thisDistance;
           matchIndex = m;
         }
+        // Fill the vector with the idxs of the flashes
+        if(thisDistance<fDistanceCandidateFlashes){
+          candidateFlashIdxs.push_back(m);
+        }
       } //End for flash
 
       //No valid match found...
@@ -461,11 +624,52 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e)
         continue;
       }
 
+      // -----------------> Here we would end with the first pool of candidate flashes <-----------------
+
+      double minChi2 = 1e6;
+      for(size_t i=0; i<candidateFlashIdxs.size(); i++)
+      {
+        fFlashTimeDebug=-99999;
+        int idx = candidateFlashIdxs[i];
+        auto & flash = flashVect[idx];
+        double _currentFlashX = flash->XCenter();
+        double _currentFlashY = flash->YCenter();
+        double _currentFlashZ = flash->ZCenter();
+        //Get the ophits associated to this opflash
+        std::vector<art::Ptr<recob::OpHit>> ophitlist = flash_ophit_assns.at(flash.key());
+        std::vector<double> ophit_z;
+        std::vector<double> ophit_y;
+        std::vector<double> ophit_weight;
+        CreateOpHitList(ophitlist, ophit_z, ophit_y, ophit_weight);
+
+        double FlashLight = this->GetFlashLight(flash->TotalPE(), direct_visibility, reflect_visibility);
+        fFlashTimeDebug = flash->Time();
+        double lightChargeRatio = (FlashLight/sliceCharge);
+        std::vector<double> LightPCA = {0.,0.};
+        GetPCA(ophit_z, ophit_y, ophit_weight, LightPCA);
+        double LightChargeCosine = abs( ChargePCA[0]*LightPCA[0] + ChargePCA[1]*LightPCA[1] ) / ( std::hypot(ChargePCA[0], ChargePCA[1]) * std::hypot(LightPCA[0], LightPCA[1]) );
+        double angle = std::acos(LightChargeCosine)*(180/M_PI);
+        double chi2;
+        if(fDo3DMatching) chi2 = std::pow(fChargeCenterX-_currentFlashX, 2)/std::pow(fXError, 2) + std::pow(fChargeCenterY-_currentFlashY, 2)/std::pow(fYError, 2) + std::pow(fChargeCenterZ-_currentFlashZ, 2)/std::pow(fZError, 2) + std::pow(angle, 2)/std::pow(fAngleError, 2);
+        else
+          chi2 = std::pow(fChargeCenterY-_currentFlashY, 2)/std::pow(fYError, 2) + std::pow(fChargeCenterZ-_currentFlashZ, 2)/std::pow(fZError, 2) + std::pow(angle, 2)/std::pow(fAngleError, 2);
+
+        if(lightChargeRatio < fLightChargeRatioBounds[0] || lightChargeRatio > fLightChargeRatioBounds[1]) continue;
+        if(chi2 < minChi2)
+        {
+          minChi2 = chi2;
+          fAngle = angle;
+          matchIndex = idx;
+        }
+      }
+
       //Best match flash pointer
       unsigned unsignedMatchIndex = matchIndex;
       const art::Ptr<recob::OpFlash> flashPtr { flashHandle, unsignedMatchIndex };
       //Update match info
       updateFlashVars(flashPtr);
+      fChi2 = minChi2;
+      fScore = 1./fChi2;
       updateMatchInfo(sliceMatchInfo);
       sliceMatchInfoVector.push_back(sliceMatchInfo);
       art::Ptr<sbn::TPCPMTBarycenterMatch> const infoPtr = makeInfoPtr(matchInfoVector->size());
@@ -476,7 +680,7 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e)
 
     int maxFlashIdx=-1;
     int minNPes = -1000000000;
-    
+
     for(int i=0; i<static_cast<int>(flashPtrVector.size()); i++)
     {
       int nPEs = flashPtrVector[i]->TotalPE();
@@ -486,6 +690,7 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e)
         minNPes = nPEs;
       }
     }
+
     if(maxFlashIdx!=-1) 
     {
       sliceAssns->addSingle(infoPtrVector[maxFlashIdx], slicePtr);
@@ -518,11 +723,12 @@ void TPCPMTBarycenterMatchProducer::InitializeSlice() {
   fFlashWidthY = -9999.;
   fFlashWidthZ = -9999.;
   fDeltaT = -9999.;
-  fDeltaX = -9999.;
   fDeltaY = -9999.;
   fDeltaZ = -9999.;
   fRadius = -9999.;
-  fDeltaX_Trigger = -9999.;
+  fChi2 = -9999.;
+  fScore = -99999.;
+  fAngle = -9999.;
   fDeltaZ_Trigger = -9999.;
   fDeltaY_Trigger = -9999.;
   fRadius_Trigger = -9999.;
@@ -538,14 +744,8 @@ void TPCPMTBarycenterMatchProducer::updateChargeVars(double sumCharge, TVector3
   fChargeWidthZ = std::sqrt( sumPosSqr[2]/sumCharge - (sumPos[2]/sumCharge)*(sumPos[2]/sumCharge) );
   fChargeTotal = sumCharge;
   if ( triggerFlashCenter[1] != -9999. ) {
-    fDeltaX_Trigger = abs(triggerFlashCenter[0] - fChargeCenterX);
-    fDeltaY_Trigger = abs(triggerFlashCenter[1] - fChargeCenterY);
-    fDeltaZ_Trigger = abs(triggerFlashCenter[2] - fChargeCenterZ);
-    if(fDo3DMatching) 
-    {
-      fRadius_Trigger = std::hypot(fDeltaX_Trigger, fDeltaY_Trigger, fDeltaZ_Trigger);
-    }
-    else
+    fDeltaY_Trigger = abs(triggerFlashCenter[0] - fChargeCenterY);
+    fDeltaZ_Trigger = abs(triggerFlashCenter[1] - fChargeCenterZ);
     fRadius_Trigger = std::hypot(fDeltaY_Trigger, fDeltaZ_Trigger);
   }
 } //End updateChargeVars()
@@ -569,12 +769,6 @@ void TPCPMTBarycenterMatchProducer::updateFlashVars(art::Ptr<recob::OpFlash> fla
   if ( fChargeT0 != -9999 ) fDeltaT = abs(matchedTime - fChargeT0);
   fDeltaY = abs(matchedYCenter - fChargeCenterY);
   fDeltaZ = abs(matchedZCenter - fChargeCenterZ);
-  if( fDo3DMatching ) 
-  {
-    fDeltaX = abs(matchedXCenter - fChargeCenterX);
-    fRadius = std::hypot(fDeltaX, fDeltaY, fDeltaZ);
-  }
-  else
   fRadius = std::hypot(fDeltaY, fDeltaZ);
 } //End updateFlashVars()
 
@@ -592,7 +786,119 @@ void TPCPMTBarycenterMatchProducer::updateMatchInfo(sbn::TPCPMTBarycenterMatch&
   matchInfo.deltaZ = fDeltaZ;
   matchInfo.radius = fRadius;
   matchInfo.radius_Trigger = fRadius_Trigger;
+  matchInfo.chi2 = fChi2;
+  matchInfo.score = fScore;
+  matchInfo.angle = fAngle;
 } //End updateMatchInfo()
 
+void TPCPMTBarycenterMatchProducer::GetPCA(std::vector<double> const&x, std::vector<double> const&y, std::vector<double> const&weight, std::vector<double> & PCA ) {
+
+  double weight_sum = std::accumulate(weight.begin(), weight.end(), 0.0);
+  //Create the matrix of points
+  size_t n = x.size();
+  Eigen::MatrixXd points(n, 2);
+  for (size_t i = 0; i < n; ++i) {
+      points(i, 0) = x[i];
+      points(i, 1) = y[i];
+  }
+
+  //Get the weighted centroid
+  Eigen::Vector2d weighted_sum(0.0, 0.0);
+  for (size_t i = 0; i < n; ++i) {
+      weighted_sum(0) += weight[i] * x[i];
+      weighted_sum(1) += weight[i] * y[i];
+  }
+
+  //Center the points wrt centroid
+  Eigen::Vector2d centroid = weighted_sum / weight_sum;
+  Eigen::MatrixXd centered_points = points.rowwise() - centroid.transpose();
+
+  //Get the weighted covariance matrix
+  Eigen::Matrix2d cov = Eigen::Matrix2d::Zero();
+  for (size_t i = 0; i < n; ++i) {
+      Eigen::Vector2d pt = centered_points.row(i);
+      cov += std::abs(weight[i]) * (pt * pt.transpose());
+  }
+  cov /= weight_sum;
+
+  // Get eigenvalues / eigenvectors 
+  Eigen::SelfAdjointEigenSolver<Eigen::Matrix2d> solver(cov);
+  Eigen::VectorXd eigvals = solver.eigenvalues();
+  Eigen::Matrix2d eigvecs = solver.eigenvectors();
+
+  // Get Principal component (largest eigenvalue)
+  int maxIndex;
+  eigvals.maxCoeff(&maxIndex);
+  Eigen::Vector2d principal_direction = eigvecs.col(maxIndex);
+
+  PCA[0] = principal_direction(0);
+  PCA[1] = principal_direction(1);
+
+}
+
+
+
+double TPCPMTBarycenterMatchProducer::GetSliceCharge(const std::vector<art::Ptr<recob::Hit>> &tpcHitsVec,  const detinfo::DetectorPropertiesData det_prop, int tpc) {
+
+    std::vector<double> plane_charge{0.,0.,0.};
+    std::vector<int>    plane_hits{0,0,0};
+    for (size_t i=0; i < tpcHitsVec.size(); i++){
+      auto hit = tpcHitsVec[i];
+      if( hit->WireID().TPC != static_cast<unsigned>(tpc) ) continue; // Skip hits not in the TPC we are considering
+      auto drift_time = (hit->PeakTime() - 500)*0.5; // assuming TPC beam readout starts at 500 ticks, conversion = 0.5 us/tick
+      double atten_correction = std::exp(drift_time/det_prop.ElectronLifetime()); // exp(us/us)
+      auto hit_plane = hit->View();
+      plane_charge.at(hit_plane) += hit->Integral()*atten_correction*(1/fCalAreaConst.at(hit_plane));
+      plane_hits.at(hit_plane)++; 
+    }
+
+    uint bestHits =  std::max_element(plane_hits.begin(), plane_hits.end()) - plane_hits.begin();
+    double _comp_charge  = plane_charge.at(bestHits);
+
+    return _comp_charge;
+}
+
+void TPCPMTBarycenterMatchProducer::CreateOpHitList( std::vector<art::Ptr<recob::OpHit>> ophitlist, std::vector<double>& ophit_z, std::vector<double>& ophit_y, std::vector<double>& ophit_weight) {
+
+  std::map<std::pair<int, int>, double> footPrintMap;
+  for (size_t i=0; i<ophitlist.size(); i++)
+  {
+    auto &ophit = ophitlist[i];
+    int channel = ophit->OpChannel();
+    int pos_z = static_cast<int>(fOpDetZ[channel]);
+    int pos_y = static_cast<int>(fOpDetY[channel]);
+    double weight = ophit->PE();
+    footPrintMap[{pos_z, pos_y}] += weight;
+  }
+
+  for (const auto& opdet : footPrintMap) {
+    int z = opdet.first.first;
+    int y = opdet.first.second;
+    int weight = opdet.second;
+    ophit_z.push_back(z);
+    ophit_y.push_back(y);
+    ophit_weight.push_back(weight*weight);
+  }
+
+}
+
+double TPCPMTBarycenterMatchProducer::GetFlashLight(double flash_pe, std::vector<double>& dir_visibility, std::vector<double>& ref_visibility) {
+
+  double tot_visibility=0;
+
+  for(size_t ch=0; ch<dir_visibility.size(); ch++){
+    if (std::find(fSkipChannelList.begin(), fSkipChannelList.end(), ch) != fSkipChannelList.end()) continue;
+    if(fOpDetType[ch]==0) tot_visibility += fOpDetVUVEff*dir_visibility[ch] + fOpDetVISEff*ref_visibility[ch];
+    else if(fOpDetType[ch]==1) tot_visibility += fOpDetVISEff*ref_visibility[ch];
+    else continue; // skip other types
+  }
+  //std::cout << " The number of PEs in the flash is " << flash_pe << " the total direct visibility is " << total_dir_visibility << " the total reflected visibility is " << total_ref_visibility << " with a VUV QE " << fOpDetVUVEff << " and a vis QE " << fOpDetVISEff << " so the total visibility is " << tot_visibility << std::endl;
+  if((flash_pe == 0) || std::isinf(1/tot_visibility))
+    return 0.0;
+  // deposited light is inverse of visibility * PE count
+  double total_gamma = (1/tot_visibility)*flash_pe;
+  return total_gamma;
+}
+
 
 DEFINE_ART_MODULE(TPCPMTBarycenterMatchProducer)
\ No newline at end of file
diff --git a/sbndcode/TPCPMTBarycenterMatching/job/sbnd_tpcpmt3dbarycentermatching_config.fcl b/sbndcode/TPCPMTBarycenterMatching/job/sbnd_tpcpmt3dbarycentermatching_config.fcl
index f6193ba3e..d6f622c8a 100644
--- a/sbndcode/TPCPMTBarycenterMatching/job/sbnd_tpcpmt3dbarycentermatching_config.fcl
+++ b/sbndcode/TPCPMTBarycenterMatching/job/sbnd_tpcpmt3dbarycentermatching_config.fcl
@@ -1,25 +1,34 @@
+#include "opticalsimparameterisations_sbnd.fcl"
+#include "digi_pmt_sbnd.fcl"
+#include "calorimetry_sbnd.fcl"
+
 BEGIN_PROLOG
 
 TPCPMTBarycenterMatchProducer:
 {
     OpFlashesModuleLabel:     ["opflashtpc0", "opflashtpc1"]
     PandoraLabel:     "pandora"
-    CollectionOnly:   true
+    CollectionOnly:   false
     Verbose:          false
     FillMatchTree:    false
     Do3DMatching:     true
+    DistanceCandidateFlashes: 100
+    VUVHits: @local::sbnd_vuv_RS100cm_hits_parameterization
+    VIVHits: @local::sbnd_vis_RS100cm_hits_parameterization
+    OpDetVUVEff: @local::sbnd_digipmt_alg.PMTCoatedVUVEff_tpc1
+    OpDetVISEff:  @local::sbnd_digipmt_alg.PMTUncoatedEff_tpc1
+    CalAreaConst: @local::sbnd_calorimetryalgdata.CalAreaConstants
+    LightChargeRatioBounds: [0.25, 3]
+    XError: 12 // cm
+    YError: 23 // cm
+    ZError: 23 // cm
+    AngleError: 27 // deg
+    FlashVetoWindow: [-1500000,1500000] // ns
     module_type: "TPCPMTBarycenterMatchProducer"
 }
 
-TPCPMTBarycenterMatchProducerSCE:
-{
-    OpFlashesModuleLabel:     ["opflashtpc0", "opflashtpc1"]
-    PandoraLabel:     "pandoraSCE"
-    CollectionOnly:   true
-    Verbose:          false
-    FillMatchTree:    false
-    Do3DMatching:     true
-    module_type: "TPCPMTBarycenterMatchProducer"
-}
+TPCPMTBarycenterMatchProducerSCE: @local::TPCPMTBarycenterMatchProducer 
+TPCPMTBarycenterMatchProducerSCE.PandoraLabel: "pandoraSCE"
+
 
 END_PROLOG
\ No newline at end of file