diff --git a/sbncode/CAFMaker/FillTrigger.cxx b/sbncode/CAFMaker/FillTrigger.cxx
index fe68e00a1..f6bb51317 100644
--- a/sbncode/CAFMaker/FillTrigger.cxx
+++ b/sbncode/CAFMaker/FillTrigger.cxx
@@ -13,6 +13,15 @@ namespace caf
     triggerInfo.global_trigger_det_time = trig.TriggerTime();
     double diff_ts = triggerInfo.global_trigger_det_time - triggerInfo.beam_gate_det_time;
     triggerInfo.trigger_within_gate = diff_ts;
+
+    triggerInfo.prev_global_trigger_time = addltrig_info.previousTriggerTimestamp;
+    triggerInfo.source_type = (int)addltrig_info.sourceType;
+    triggerInfo.trigger_type = (int)addltrig_info.triggerType;
+    triggerInfo.trigger_id = addltrig_info.triggerID;
+    triggerInfo.gate_id = addltrig_info.gateID;
+    triggerInfo.trigger_count = addltrig_info.triggerCount;
+    triggerInfo.gate_count = addltrig_info.gateCount;
+    triggerInfo.gate_delta = (int)addltrig_info.gateCountFromPreviousTrigger;
   }
 
   void FillTriggerMC(double absolute_time, caf::SRTrigger& triggerInfo) {
diff --git a/sbncode/CAFMaker/FillTrue.cxx b/sbncode/CAFMaker/FillTrue.cxx
index 3b9da6340..5d7393832 100644
--- a/sbncode/CAFMaker/FillTrue.cxx
+++ b/sbncode/CAFMaker/FillTrue.cxx
@@ -470,7 +470,9 @@ namespace caf {
     for(const caf::SRTrueParticle& part: srparticles){
       // save the G4 particles that came from this interaction
       if(part.interaction_id == (int)i) {
-        if(part.start_process == caf::kG4primary) srneutrino.prim.push_back(part);
+        //if(part.start_process == caf::kG4primary) srneutrino.prim.push_back(part);
+        // Feb. 15th 2024, Jaesung Kim: Saving non-primaries for G4Reweight study
+        srneutrino.prim.push_back(part);
 
         // total up the deposited energy
         for(int p = 0; p < 3; ++p) { 
diff --git a/sbncode/Calibration/TrackCaloSkimmer_module.cc b/sbncode/Calibration/TrackCaloSkimmer_module.cc
index 8cbfcaa19..b31b86c9e 100644
--- a/sbncode/Calibration/TrackCaloSkimmer_module.cc
+++ b/sbncode/Calibration/TrackCaloSkimmer_module.cc
@@ -831,6 +831,41 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle,
     }
   }
 
+  // Save the true trajectory
+  for (unsigned i_traj = 0; i_traj < particle.NumberTrajectoryPoints(); i_traj++) {
+    // Get trajectory point
+    TVector3 traj = particle.Position(i_traj).Vect();
+    geo::Point_t traj_p(traj.X(), traj.Y(), traj.Z());
+
+    // lookup TPC
+    geo::TPCID tpc; // invalid by default
+    for (auto const &cryo: geo->Iterate<geo::CryostatGeo>()) {
+      for (auto const& TPC : geo->Iterate<geo::TPCGeo>(cryo.ID())) {
+        if (TPC.ActiveBoundingBox().ContainsPosition(traj_p)) {
+          tpc = TPC.ID();
+          break;
+        }
+      }
+      if (tpc.isValid) break;
+    }
+
+    // add in space-charge-deflected position if applicable
+    geo::Point_t traj_p_sce = tpc.isValid ? TrajectoryToWirePosition(traj_p, tpc) : traj_p;
+
+    sbn::Vector3D traj_v;
+    traj_v.x = traj_p.x();
+    traj_v.y = traj_p.y();
+    traj_v.z = traj_p.z();
+
+    sbn::Vector3D traj_v_sce;
+    traj_v_sce.x = traj_p_sce.x();
+    traj_v_sce.y = traj_p_sce.y();
+    traj_v_sce.z = traj_p_sce.z();
+
+    trueparticle.traj.push_back(traj_v);
+    trueparticle.traj_sce.push_back(traj_v_sce);
+  }
+
   return trueparticle;
 }
 
diff --git a/sbncode/EventGenerator/MeVPrtl/Tools/ALP/Meson2ALP_tool.cc b/sbncode/EventGenerator/MeVPrtl/Tools/ALP/Meson2ALP_tool.cc
index fad80f2a1..2f709e026 100644
--- a/sbncode/EventGenerator/MeVPrtl/Tools/ALP/Meson2ALP_tool.cc
+++ b/sbncode/EventGenerator/MeVPrtl/Tools/ALP/Meson2ALP_tool.cc
@@ -76,6 +76,7 @@ class Meson2ALP : public IMeVPrtlFlux
   double fcG; //!< Axion coupling to gluon
   double fcB; //!< Axion coupling to U(1) B boson (before EW symmetry breaking)
   double fcW; //!< Axion coupling to SU(2) W bosons (before EW symmetry breaking)
+  bool fIncludeMassSuppressionFactor; // Whether to include a phenomonological suppression to mixing at high mass
 
   // branching ratios
   double fPi0BR;
@@ -115,6 +116,7 @@ void Meson2ALP::configure(fhicl::ParameterSet const &pset)
   fMaxWeightPi0 = pset.get<double>("MaxWeightPi0", 1.);
   fMaxWeightEta = pset.get<double>("MaxWeightEta", 1.);
   fMaxWeightEtaP = pset.get<double>("MaxWeightEtaP", 1.);
+  fIncludeMassSuppressionFactor = pset.get<bool>("IncludeMassSuppressionFactor", false);
 
   fPi0BR = Pi0BR();
   fEtaBR = EtaBR();
@@ -136,7 +138,7 @@ double Meson2ALP::EtaBR() const {
 
   double mixing_angle = (1. / sqrt(6)) * (fcG * fpion / ffa) * (fM*fM - (4./9.)*pzero_mass*pzero_mass) / (fM*fM - eta_mass*eta_mass);
   double qcd_rate_f = 1.;
-  if (fM > eta_mass) qcd_rate_f = pow(fM/eta_mass, -1.6);
+  if (fIncludeMassSuppressionFactor && fM > eta_mass) qcd_rate_f = pow(fM/eta_mass, -1.6);
 
   return mixing_angle*mixing_angle*qcd_rate_f;
 }
@@ -148,7 +150,7 @@ double Meson2ALP::EtaPBR() const {
 
   double mixing_angle = (1. / sqrt(12)) * (fcG * fpion / ffa) * (fM*fM - (16./9.)*pzero_mass*pzero_mass) / (fM*fM - etap_mass*etap_mass);
   double qcd_rate_f = 1.;
-  if (fM > etap_mass) qcd_rate_f = pow(fM/etap_mass, -1.6);
+  if (fIncludeMassSuppressionFactor && fM > etap_mass) qcd_rate_f = pow(fM/etap_mass, -1.6);
 
   return mixing_angle*mixing_angle*qcd_rate_f;
 }
@@ -160,7 +162,7 @@ double Meson2ALP::Pi0BR() const {
   double mixing_angle = (1./6) * (fcG * fpion / ffa) * fM*fM / (fM*fM - pzero_mass*pzero_mass);
 
   double qcd_rate_f = 1.;
-  if (fM > pzero_mass) qcd_rate_f = pow(fM/pzero_mass, -1.6);
+  if (fIncludeMassSuppressionFactor && fM > pzero_mass) qcd_rate_f = pow(fM/pzero_mass, -1.6);
 
   return mixing_angle*mixing_angle*qcd_rate_f;
 }
@@ -177,6 +179,10 @@ bool Meson2ALP::MakeFlux(const simb::MCFlux &flux, evgen::ldm::MeVPrtlFlux &alp,
 
   // Energy is same as for meson (don't worry about momentum conservation)
   double alp_energy = meson.mom.E();
+
+  // ignore if alp isn't energetic enough
+  if (alp_energy < fM) return false;
+
   double alp_momentum = sqrt(alp_energy*alp_energy - fM*fM);
 
   // Momentum 4-vector
diff --git a/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx b/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx
index 8ae583973..caf10e17f 100644
--- a/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx
+++ b/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx
@@ -1,6 +1,9 @@
 #include "TrajectoryMCSFitter.h"
 #include "lardataobj/RecoBase/Track.h"
 #include "larcorealg/Geometry/geo_vectors_utils.h"
+#include "larcore/Geometry/Geometry.h"
+#include "larcore/CoreUtils/ServiceUtil.h"
+#include "larcorealg/Geometry/GeometryCore.h"
 #include "TMatrixDSym.h"
 #include "TMatrixDSymEigen.h"
 
@@ -15,7 +18,8 @@ recob::MCSFitResult TrajectoryMCSFitter::fitMcs(const recob::TrackTrajectory& tr
   vector<size_t> breakpoints;
   vector<float> segradlengths;
   vector<float> cumseglens;
-  breakTrajInSegments(traj, breakpoints, segradlengths, cumseglens);
+  vector<bool> breakpointsgood;
+  breakTrajInSegments(traj, breakpoints, segradlengths, cumseglens, breakpointsgood);
   //
   // Fit segment directions, and get 3D angles between them
   //
@@ -28,7 +32,11 @@ recob::MCSFitResult TrajectoryMCSFitter::fitMcs(const recob::TrackTrajectory& tr
     if (p>0) {
       if (segradlengths[p]<-100. || segradlengths[p-1]<-100.) {
         dtheta.push_back(-999.);
-            } else { 
+      } 
+      if (!breakpointsgood[p] || !breakpointsgood[p-1]) {
+        dtheta.push_back(-999.);
+      }
+      else { 
         const double cosval = pcdir0.X()*pcdir1.X()+pcdir0.Y()*pcdir1.Y()+pcdir0.Z()*pcdir1.Z();
         //assert(std::abs(cosval)<=1);
         //units are mrad
@@ -49,15 +57,21 @@ recob::MCSFitResult TrajectoryMCSFitter::fitMcs(const recob::TrackTrajectory& tr
   }
   const ScanResult fwdResult = doLikelihoodScan(dtheta, segradlengths, cumLenFwd, true,  momDepConst, pid);
   const ScanResult bwdResult = doLikelihoodScan(dtheta, segradlengths, cumLenBwd, false, momDepConst, pid);
-  //
+  // std::cout << "fwdResult.p=" << fwdResult.p << " fwdResult.pUnc=" << fwdResult.pUnc << " fwdResult.logL=" << fwdResult.logL << std::endl;
   return recob::MCSFitResult(pid,
                             fwdResult.p,fwdResult.pUnc,fwdResult.logL,
                             bwdResult.p,bwdResult.pUnc,bwdResult.logL,
                             segradlengths,dtheta);
 }
 
-void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj, vector<size_t>& breakpoints, vector<float>& segradlengths, vector<float>& cumseglens) const {
-  //
+void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj, vector<size_t>& breakpoints, vector<float>& segradlengths, vector<float>& cumseglens, vector<bool>& breakpointsgood) const {
+  // set fiducial volume
+  std::vector<geo::BoxBoundedGeo> fiducialVolumes;
+  fiducialVolumes = setFiducialVolumes();
+  // set excluded volumes
+  std::vector<geo::BoxBoundedGeo> excludeVolumes;
+  excludeVolumes = setExcludeVolumes();
+
   const double trajlen = traj.Length();
   const int nseg = std::max(minNSegs_,int(trajlen/segLen_));
   const double thisSegLen = trajlen/double(nseg);
@@ -69,6 +83,8 @@ void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj
   auto nextValid=traj.FirstValidPoint();
   breakpoints.push_back(nextValid);
   auto pos0 = traj.LocationAtPoint(nextValid);
+  bool pos0good = isInVolume(fiducialVolumes, pos0) && !isInVolume(excludeVolumes, pos0);
+  breakpointsgood.push_back(pos0good);
   nextValid = traj.NextValidPoint(nextValid+1);
   int npoints = 0;
   while (nextValid!=recob::TrackTrajectory::InvalidIndex) {
@@ -77,7 +93,9 @@ void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj
     pos0=pos1;
     npoints++;
     if (thislen>=thisSegLen) {
+      bool pos1good = isInVolume(fiducialVolumes, pos1) && !isInVolume(excludeVolumes, pos1);
       breakpoints.push_back(nextValid);
+      breakpointsgood.push_back(pos1good);
       if (npoints>=minHitsPerSegment_) segradlengths.push_back(thislen*lar_radl_inv);
       else segradlengths.push_back(-999.);
       cumseglens.push_back(cumseglens.back()+thislen);
@@ -88,7 +106,10 @@ void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj
   }
   //then add last segment
   if (thislen>0.) {
+    auto endpointpos = traj.LocationAtPoint(traj.LastValidPoint());
+    bool endpointposgood = isInVolume(fiducialVolumes, endpointpos) && !isInVolume(excludeVolumes, endpointpos);
     breakpoints.push_back(traj.LastValidPoint()+1);
+    breakpointsgood.push_back(endpointposgood);
     segradlengths.push_back(thislen*lar_radl_inv);
     cumseglens.push_back(cumseglens.back()+thislen);
   }
@@ -202,7 +223,6 @@ double TrajectoryMCSFitter::mcsLikelihood(double p, double theta0x, std::vector<
   double result = 0;
   for (int i = beg; i != end; i+=incr ) {
     if (dthetaij[i]<0) {
-      //cout << "skip segment with too few points" << endl;
       continue;
     }
     //
@@ -302,3 +322,75 @@ double TrajectoryMCSFitter::GetE(const double initial_E, const double length_tra
   }
   return current_E;
 }
+
+std::vector<geo::BoxBoundedGeo> TrajectoryMCSFitter::setFiducialVolumes() const {
+  std::vector<geo::BoxBoundedGeo> fiducialVolumes;
+  std::vector<std::vector<geo::BoxBoundedGeo>> TPCVolumes;
+  std::vector<geo::BoxBoundedGeo> ActiveVolumes;
+  
+  const geo::GeometryCore *geometry = lar::providerFrom<geo::Geometry>();
+  for (auto const &cryoID: geometry->Iterate<geo::CryostatID>()) {
+    std::vector<geo::BoxBoundedGeo> thisTPCVolumes;
+    for (auto const& tpc: geometry->Iterate<geo::TPCGeo>(cryoID)) {
+      thisTPCVolumes.push_back(tpc.ActiveBoundingBox());
+    }
+    TPCVolumes.push_back(std::move(thisTPCVolumes));
+  }
+  for (const std::vector<geo::BoxBoundedGeo> &tpcs: TPCVolumes) {
+    double xMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinX() < rhs.MinX(); })->MinX();
+    double xMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxX() < rhs.MaxX(); })->MaxX();
+    double yMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinY() < rhs.MinY(); })->MinY();
+    double yMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxY() < rhs.MaxY(); })->MaxY();
+    double zMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinZ() < rhs.MinZ(); })->MinZ();
+    double zMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxZ() < rhs.MaxZ(); })->MaxZ();
+    ActiveVolumes.emplace_back(xMin, xMax, yMin, yMax, zMin, zMax);
+  }
+  double fidInsetMinX = 0;
+  double fidInsetMaxX = 0;
+  double fidInsetMinY = 0;
+  double fidInsetMaxY = 0;
+  double fidInsetMinZ = 0;
+  double fidInsetMaxZ = 0;
+
+  if (fiducialVolumeInsets_.size() == 6) {
+    fidInsetMinX = fiducialVolumeInsets_[0];
+    fidInsetMaxX = fiducialVolumeInsets_[1];
+    fidInsetMinY = fiducialVolumeInsets_[2];
+    fidInsetMaxY = fiducialVolumeInsets_[3];
+    fidInsetMinZ = fiducialVolumeInsets_[4];
+    fidInsetMaxZ = fiducialVolumeInsets_[5];
+  }
+  else if (!fiducialVolumeInsets_.empty()) {
+    std::cout << "Error: fiducialVolumeInsets vector must have length of 6, not fiducializing" << std::endl;
+  }
+  for (const geo::BoxBoundedGeo &AV: ActiveVolumes) {
+    fiducialVolumes.emplace_back(AV.MinX() + fidInsetMinX, AV.MaxX() - fidInsetMaxX,
+                                 AV.MinY() + fidInsetMinY, AV.MaxY() - fidInsetMaxY,
+                                 AV.MinZ() + fidInsetMinZ, AV.MaxZ() - fidInsetMaxZ);
+  }
+  return fiducialVolumes;
+}
+
+std::vector<geo::BoxBoundedGeo> TrajectoryMCSFitter::setExcludeVolumes() const {
+  std::vector<geo::BoxBoundedGeo> excludeVolumes;
+  if (excludeVolumes_.size()%6 != 0) {
+    std::cout << "Error: excludeVolumes vector must have length multiple of 6, not excluding any regions" << std::endl;
+    excludeVolumes.emplace_back(-9999.0, -9999.0, -9999.0, -9999.0, -9999.0, -9999.0);
+    return excludeVolumes;
+  }
+  for (unsigned int i=0; i<excludeVolumes_.size()/6; i++) {
+    excludeVolumes.emplace_back(excludeVolumes_[6*i+0], excludeVolumes_[6*i+1], excludeVolumes_[6*i+2], excludeVolumes_[6*i+3], excludeVolumes_[6*i+4], excludeVolumes_[6*i+5]);
+  }
+  return excludeVolumes;
+}
+
+bool TrajectoryMCSFitter::isInVolume(const std::vector<geo::BoxBoundedGeo> &volumes, const geo::Point_t &point) const {
+  for (const geo::BoxBoundedGeo &volume: volumes) {
+    if (point.X()>=volume.MinX() && point.X()<=volume.MaxX() &&
+        point.Y()>=volume.MinY() && point.Y()<=volume.MaxY() &&
+        point.Z()>=volume.MinZ() && point.Z()<=volume.MaxZ()) {
+      return true;
+    }
+  }
+  return false;
+}
diff --git a/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h b/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h
index 07b63c443..6d3b48032 100644
--- a/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h
+++ b/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h
@@ -3,11 +3,15 @@
 
 #include "fhiclcpp/ParameterSet.h"
 #include "fhiclcpp/types/Atom.h"
+#include "fhiclcpp/types/Sequence.h"
 #include "fhiclcpp/types/Table.h"
 #include "canvas/Persistency/Common/Ptr.h"
 #include "lardataobj/RecoBase/MCSFitResult.h"
 #include "lardataobj/RecoBase/Track.h"
 #include "lardata/RecoObjects/TrackState.h"
+#include "larcore/Geometry/Geometry.h"
+#include "larcore/CoreUtils/ServiceUtil.h"
+#include "larcorealg/Geometry/GeometryCore.h"
 
 namespace trkf::sbn {
   /**
@@ -87,10 +91,18 @@ namespace trkf::sbn {
         Comment("Angular resolution parameter used in modified Highland formula. Unit is mrad."),
         3.0
       };
+      fhicl::Sequence<double> fiducialVolumeInsets {
+        Name("fiducialVolumeInsets"),
+        Comment("insets from the TPC boundaries to exclude from the fit")
+      };
+      fhicl::Sequence<double> excludeVolumes {
+        Name("excludeVolumes"),
+        Comment("regions to exclude")
+      };
     };
     using Parameters = fhicl::Table<Config>;
     //
-    TrajectoryMCSFitter(int pIdHyp, int minNSegs, double segLen, int minHitsPerSegment, int nElossSteps, int eLossMode, double pMin, double pMax, double pStep, double angResol){
+    TrajectoryMCSFitter(int pIdHyp, int minNSegs, double segLen, int minHitsPerSegment, int nElossSteps, int eLossMode, double pMin, double pMax, double pStep, double angResol, std::vector<double>fiducialVolumeInsets, std::vector<double> excludeVolumes){
       pIdHyp_ = pIdHyp;
       minNSegs_ = minNSegs;
       segLen_ = segLen;
@@ -101,9 +113,11 @@ namespace trkf::sbn {
       pMax_ = pMax;
       pStep_ = pStep;
       angResol_ = angResol;
+      fiducialVolumeInsets_ = fiducialVolumeInsets;
+      excludeVolumes_ = excludeVolumes;
     }
     explicit TrajectoryMCSFitter(const Parameters & p)
-      : TrajectoryMCSFitter(p().pIdHypothesis(),p().minNumSegments(),p().segmentLength(),p().minHitsPerSegment(),p().nElossSteps(),p().eLossMode(),p().pMin(),p().pMax(),p().pStep(),p().angResol()) {}
+      : TrajectoryMCSFitter(p().pIdHypothesis(),p().minNumSegments(),p().segmentLength(),p().minHitsPerSegment(),p().nElossSteps(),p().eLossMode(),p().pMin(),p().pMax(),p().pStep(),p().angResol(),p().fiducialVolumeInsets(),p().excludeVolumes()) {}
     //
     recob::MCSFitResult fitMcs(const recob::TrackTrajectory& traj, bool momDepConst = true) const { return fitMcs(traj,pIdHyp_,momDepConst); }
     recob::MCSFitResult fitMcs(const recob::Track& track,          bool momDepConst = true) const { return fitMcs(track,pIdHyp_,momDepConst); }
@@ -117,7 +131,7 @@ namespace trkf::sbn {
       return fitMcs(tt,pid,momDepConst);
     }
     //
-    void breakTrajInSegments(const recob::TrackTrajectory& traj, std::vector<size_t>& breakpoints, std::vector<float>& segradlengths, std::vector<float>& cumseglens) const;
+    void breakTrajInSegments(const recob::TrackTrajectory& traj, std::vector<size_t>& breakpoints, std::vector<float>& segradlengths, std::vector<float>& cumseglens, std::vector<bool>& breakpointsgood) const;
     void linearRegression(const recob::TrackTrajectory& traj, const size_t firstPoint, const size_t lastPoint, recob::tracking::Vector_t& pcdir) const;
     double mcsLikelihood(double p, double theta0x, std::vector<float>& dthetaij, std::vector<float>& seg_nradl, std::vector<float>& cumLen, bool fwd, bool momDepConst, int pid) const;
     //
@@ -146,6 +160,9 @@ namespace trkf::sbn {
     double energyLossLandau(const double mass2,const double E2, const double x) const;
     //
     double GetE(const double initial_E, const double length_travelled, const double mass) const;
+    std::vector<geo::BoxBoundedGeo> setFiducialVolumes() const;
+    std::vector<geo::BoxBoundedGeo> setExcludeVolumes() const;
+    bool isInVolume(const std::vector<geo::BoxBoundedGeo> &volumes, const geo::Point_t &point) const;
     //
   private:
     int    pIdHyp_;
@@ -158,6 +175,8 @@ namespace trkf::sbn {
     double pMax_;
     double pStep_;
     double angResol_;
+    std::vector<double> fiducialVolumeInsets_;
+    std::vector<double> excludeVolumes_;
   };
 }
 
diff --git a/sbncode/LArRecoProducer/mcsproducer.fcl b/sbncode/LArRecoProducer/mcsproducer.fcl
index e41a6d00b..5db98243c 100644
--- a/sbncode/LArRecoProducer/mcsproducer.fcl
+++ b/sbncode/LArRecoProducer/mcsproducer.fcl
@@ -1,7 +1,10 @@
 BEGIN_PROLOG
 mcs_sbn: {
   module_type: MCSFitAllPID
-  MCS: {}
+  MCS: {
+    fiducialVolumeInsets: [] 
+    excludeVolumes: [] 
+  }
   TrackLabel: pandoraTrack
 }
 
diff --git a/sbncode/LArRecoProducer/mcsproducer_icarus.fcl b/sbncode/LArRecoProducer/mcsproducer_icarus.fcl
new file mode 100644
index 000000000..4926fdfc7
--- /dev/null
+++ b/sbncode/LArRecoProducer/mcsproducer_icarus.fcl
@@ -0,0 +1,13 @@
+BEGIN_PROLOG
+mcs_sbn: {
+  module_type: MCSFitAllPID
+  MCS: {
+    eLossMode: 2
+    angResol: 10.0
+    fiducialVolumeInsets: [10.0, 10.0, 10.0, 10.0, 10.0, 10.0] 
+    excludeVolumes: [190.0, 240.0, -9999.0, 9999.0, -9999.0, 9999.0] 
+  }
+  TrackLabel: pandoraTrack
+}
+
+END_PROLOG
diff --git a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc
index 4793171b9..6e836054c 100644
--- a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc
+++ b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc
@@ -116,19 +116,22 @@ void SystToolsEventWeight::produce(art::Event& e) {
           << sp->GetFullyQualifiedName() << "\n";
       }
 
-      //==== syst_resp->size() is (Number of MCTruth)
-      //==== Each index corresponds to each of MCTruth
+      // syst_resp->size() is (Number of MCTruth)
+      // Each index corresponds to each of MCTruth
       int nMCTruthIndex = 0;
       if(fDebugMode) std::cout << "[SystToolsEventWeight::produce]   syst_resp.size() (= Number of MCTruth) of this SystProvider = " << syst_resp->size() << std::endl;
 
-      //==== Looping over syst_resp is identical to looping over MCTruth
+      // Looping over syst_resp is identical to looping over MCTruth
       for(systtools::event_unit_response_t const& resp: *syst_resp) {
-        //==== resp.size() corresponds to number of knobs we altered;
-        //==== e.g., MaCCQE, MaCCRES, MvCCRE -> resp.size() = 3
+        // resp.size() corresponds to number of knobs we altered;
+        // e.g., MaCCQE, MaCCRES, MvCCRE -> resp.size() = 3
         if(fDebugMode){
           std::cout << "[SystToolsEventWeight::produce]     sp->GetSystMetaData().size() (expected) = " << sp->GetSystMetaData().size() << "\n"
                     << "[SystToolsEventWeight::produce]     resp.size() of this syst_resp (produced) = " << resp.size() << std::endl;
         }
+
+        // Below check is not valid if there is a dependent dial
+        /*
         if(sp->GetSystMetaData().size()!=resp.size()){
           throw cet::exception{ "SystToolsEventWeight" } 
             << "sp->GetFullyQualifiedName() = " << sp->GetFullyQualifiedName() << "\n"
@@ -136,22 +139,32 @@ void SystToolsEventWeight::produce(art::Event& e) {
                     << resp.size() << " are produced. "
                     << "Probably this particular event is not relevant to this systematic variation.\n";
         }
+        */
 
         sbn::evwgh::EventWeightMap mcwgh;
 
         for( auto const& r: resp){
 
-          //==== responses.size() is the number of universes
+          // responses.size() is the number of universes
           systtools::SystParamHeader const &sph = fParamHeaderHelper.GetHeader( r.pid );
           std::string prettyName = sph.prettyName;
 
+          if(sph.isResponselessParam) continue;
+
           if(fDebugMode){
             std::cout << "[SystToolsEventWeight::produce]       pid of this resp = " << r.pid << "\n"
                       << "[SystToolsEventWeight::produce]       prettyName of this resp = " << prettyName << "\n"
                       << "[SystToolsEventWeight::produce]       paramVariations.size() of this resp (expected) = " << sph.paramVariations.size() << "\n"
                       << "[SystToolsEventWeight::produce]       responses.size() of this resp (produced) = " << r.responses.size() << std::endl;
           }
-          if(sph.paramVariations.size()!=r.responses.size()){
+
+          if(sph.isCorrection && (r.responses.size()!=1)){
+            throw cet::exception{ "SystToolsEventWeight" }
+              << "prettyName of this resp = " << prettyName << "\n"
+              << "We expect to have 1 universes as it is a correction, but "
+                      << r.responses.size() << " are produced\n";
+          }
+          if(!sph.isCorrection && sph.paramVariations.size()!=r.responses.size()){
             throw cet::exception{ "SystToolsEventWeight" }
               << "prettyName of this resp = " << prettyName << "\n"
               << "We expect to have " << sph.paramVariations.size() << " universes, but "
@@ -159,8 +172,8 @@ void SystToolsEventWeight::produce(art::Event& e) {
                       << "Probably this particular event is not relevant to this systematic variation.\n";
           }
 
-          //==== r.responses : std::vector<double>
-          //==== std::map<std::string, std::vector<float> > EventWeightMap
+          // r.responses : std::vector<double>
+          // std::map<std::string, std::vector<float> > EventWeightMap
           mcwgh[sp->GetFullyQualifiedName()+"_"+prettyName].assign(r.responses.cbegin(), r.responses.cend());
 
         } // END loop over knobs
@@ -191,25 +204,92 @@ void SystToolsEventWeight::beginRun(art::Run& run) {
                                   << "sp->GetFullyQualifiedName() = " << sp->GetFullyQualifiedName() << "\n"
                                   << "sp->GetInstanceName() = " << sp->GetInstanceName() << "\n"
                                   << "Printing each SystParamHeader of this ISystProviderTool..";
-    //==== Note: typedef std::vector<SystParamHeader> SystMetaData;
+    // Note: typedef std::vector<SystParamHeader> SystMetaData;
     auto const& smd = sp->GetSystMetaData();
+
+    // make a map of responsless-response params
+    std::map<systtools::paramId_t, std::vector<systtools::paramId_t>> map_resp_to_respless;
     for( auto &sph : smd ){
-      if (fDebugMode) {
-        std::cout << "  sph.prettyName = " << sph.prettyName << std::endl;
+      if(sph.isResponselessParam){
+        auto it = map_resp_to_respless.find( sph.responseParamId );
+        if( it != map_resp_to_respless.end() ){
+          it->second.push_back( sph.systParamId );
+        }
+        else{
+          map_resp_to_respless[sph.responseParamId] = {};
+          map_resp_to_respless[sph.responseParamId].push_back( sph.systParamId );
+        }
       }
+    }
+    if (fDebugMode) {
+      for(const auto& it: map_resp_to_respless){
+        const auto& sph = systtools::GetParam(smd, it.first);
+        std::cout << "Found a dependent dial: " << sph.prettyName << std::endl;
+        for(const auto& depdialid: it.second){
+          const auto& sph_dep = systtools::GetParam(smd, depdialid);
+          std::cout << "- dep dial: " << sph_dep.prettyName << std::endl;
+        }
+      }
+    }
 
-      std::string rwmode = "multisigma";
-      if(sph.isRandomlyThrown) rwmode = "multisim";
+    for( auto &sph : smd ){
 
-      if (fDebugMode) {
-        std::cout << "  rwmode = " << rwmode << std::endl;
+      // responsless
+      if(sph.isResponselessParam){
+        if (fDebugMode) {
+          std::cout << "Responsless dial found: " << sph.prettyName << ", thus skipping" << std::endl;
+        }
+        continue;
       }
 
-      std::vector<float> widths { sph.paramVariations.begin(), sph.paramVariations.end() };
-
       sbn::evwgh::EventWeightParameterSet fParameterSet;
-      fParameterSet.AddParameter(sph.prettyName, std::move(widths));
-      fParameterSet.Configure(sp->GetFullyQualifiedName()+"_"+sph.prettyName, rwmode, sph.paramVariations.size());
+      std::string rwmode = "";
+
+      std::vector<double> paramVars = sph.isCorrection ? std::vector<double>(1, sph.centralParamValue) : sph.paramVariations;
+
+      auto it = map_resp_to_respless.find( sph.systParamId );
+      if(it!=map_resp_to_respless.end()){
+
+        for(const auto depdialid: it->second){
+          const auto& sph_dep = systtools::GetParam(smd, depdialid);
+          std::vector<double> paramVars_dep = sph_dep.isCorrection ? std::vector<double>(1, sph_dep.centralParamValue) : sph_dep.paramVariations;
+
+          if(rwmode=="") rwmode = sph_dep.isRandomlyThrown ? "multisim" : "multisigma";
+          else{
+            if(rwmode!= (sph_dep.isRandomlyThrown ? "multisim" : "multisigma")){
+              throw cet::exception{ "SystToolsEventWeight" }
+                << sph.prettyName << " depends on other dials, but the remode are different between the deps dials\n";
+            }
+          }
+
+          std::vector<float> widths_dep { paramVars_dep.begin(), paramVars_dep.end() };
+          fParameterSet.AddParameter(sph_dep.prettyName, std::move(widths_dep));
+
+        }
+
+
+
+      }
+      else{
+
+        rwmode = sph.isRandomlyThrown ? "multisim" : "multisigma";
+
+        std::vector<float> widths { paramVars.begin(), paramVars.end() };
+
+        fParameterSet.AddParameter(sph.prettyName, std::move(widths));
+
+      }
+
+      if(rwmode==""){
+        throw cet::exception{ "SystToolsEventWeight" }
+          << "rwmode not set for " << sph.prettyName << "\n";
+      }
+
+      if (fDebugMode) {
+        std::cout << "  sph.prettyName = " << sph.prettyName << ", rwmode = " << rwmode << " is added to the header" << std::endl;
+      }
+
+      fParameterSet.Configure(sp->GetFullyQualifiedName()+"_"+sph.prettyName, rwmode, paramVars.size());
       fParameterSet.FillKnobValues();
 
       p->push_back( std::move(fParameterSet) );
diff --git a/sbncode/TPCReco/VertexStub/config/sbn_stub_merge_tools.fcl b/sbncode/TPCReco/VertexStub/config/sbn_stub_merge_tools.fcl
index ec4bb6a13..78e14279f 100644
--- a/sbncode/TPCReco/VertexStub/config/sbn_stub_merge_tools.fcl
+++ b/sbncode/TPCReco/VertexStub/config/sbn_stub_merge_tools.fcl
@@ -13,6 +13,8 @@ twoplane_stub_merge: {
   SaveOldStubs: false
 }
 
-stub_merge: [@local::plane_stub_merge, @local::twoplane_stub_merge]
+# stub_merge: [@local::plane_stub_merge, @local::twoplane_stub_merge]
+# NOTE: disable multi-plane matching while other planes are not calibrated
+stub_merge: [@local::plane_stub_merge]
 
 END_PROLOG
diff --git a/ups/product_deps b/ups/product_deps
index 8894b965c..fcd192934 100644
--- a/ups/product_deps
+++ b/ups/product_deps
@@ -257,10 +257,10 @@ larcv2			v2_1_0		-
 larsoft			v09_72_00_01	-
 sbnanaobj		v09_20_06_03	-
 sbndaq_artdaq_core	v1_06_00of0	-
-sbndata			v01_04		-
+sbndata			v01_05		-
 sbnobj			v09_16_00_02	-
-systematicstools	v01_02_00	-
-nusystematics		v01_02_10	-
+systematicstools	v01_02_01	-
+nusystematics		v01_02_11	-
 cetmodules		v3_20_00	-	only_for_build
 end_product_list
 ####################################