From 8265839752040f071149dbb565bbde1790489563 Mon Sep 17 00:00:00 2001
From: michaelmackenzie <michaelmackenzie@u.northwestern.edu>
Date: Tue, 13 Oct 2026 14:44:35 -0500
Subject: [PATCH 1/5] HelixTool: don't count straws without signal in the helix
 hit ratio

The hit ratio (measured straw hits / expected plane-face crossings) counted
every face the helix crosses inside the tracker active radius as an
expected hit, so dead panels and straws (TrackerStatus noSignal/suppress)
biased it low and made minHitRatio cuts depend on the tracker conditions.

- HelixTool takes an optional StrawNoSignal function. With it, each face
  crossing is weighted by the live fraction of the nearest straw in each
  layer at the crossing point. A function is used instead of TrackerStatus
  to avoid a Mu2eUtilities -> TrackerConditions library dependency.
- The hit ratio is evaluated lazily, on the first hitRatio() call; the
  filters apply it as the last cut, and skip it when minHitRatio <= 0.
- HelixFilter and MultiHelixFilter: new useTrackerStatus parameter
  (default true) passes TrackerStatus noSignal || suppress; false keeps
  the previous behavior. ReadTriggerInfo uses the same definition.
- Removed the unused HelixTool::nMinHitsLoop() (never set) and
  _trackerLength; nstrawhits() returns int.

With the Sim_commissioning reduced-tracker conditions (31% of straws
without signal) the median CE helix hit ratio goes from 0.56 to 1.0, so
minHitRatio 0.4 keeps its meaning. With useTrackerStatus false the
trigger results are identical to before.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 Mu2eUtilities/inc/HelixTool.hh            | 31 +++++---
 Mu2eUtilities/src/HelixTool.cc            | 94 +++++++++++++++++++----
 Trigger/CMakeLists.txt                    |  2 +
 Trigger/src/ReadTriggerInfo_module.cc     |  9 ++-
 Trigger/src/SConscript                    |  1 +
 TrkFilters/CMakeLists.txt                 |  4 +
 TrkFilters/src/HelixFilter_module.cc      | 27 +++++--
 TrkFilters/src/MultiHelixFilter_module.cc | 25 ++++--
 TrkFilters/src/SConscript                 |  2 +
 9 files changed, 158 insertions(+), 37 deletions(-)

diff --git a/Mu2eUtilities/inc/HelixTool.hh b/Mu2eUtilities/inc/HelixTool.hh
index 7c4c19e567..69582ab65d 100644
--- a/Mu2eUtilities/inc/HelixTool.hh
+++ b/Mu2eUtilities/inc/HelixTool.hh
@@ -4,6 +4,8 @@
 // Original author G. Pezzullo
 //
 
+#include <functional>
+
 #include "Offline/DataProducts/inc/StrawId.hh"
 #include "Math/GenVector/VectorUtil.h"
 #include "Math/VectorUtil.h"
@@ -13,6 +15,7 @@
 #include "Offline/TrackerGeom/inc/Tracker.hh"
 
 namespace mu2e { class Tracker; }
+namespace mu2e { class Plane; }
 namespace mu2e { struct HelixSeed; }
 
 using namespace ROOT::Math::VectorUtil;
@@ -22,35 +25,43 @@ namespace mu2e {
   class HelixTool{
 
   public:
-    HelixTool(const HelixSeed *Helix, const mu2e::Tracker*MyTracker);
+    // returns true for a straw that can't produce a hit (e.g. TrackerStatus noSignal or suppress).
+    // An empty function means all straws are live. A function is used instead of TrackerStatus
+    // to avoid a library dependency of Mu2eUtilities on TrackerConditions.
+    using StrawNoSignal = std::function<bool(StrawId const&)>;
+
+    HelixTool(const HelixSeed *Helix, const mu2e::Tracker*MyTracker, StrawNoSignal NoSignal = StrawNoSignal());
 
     // Accept compiler supplied d'tor, copy c'tor and assignment operator.
 
     // Accessors
     float  nLoops           () const { return _nLoops;            }
-    int    nMinHitsLoop     () const { return _nMinHitsLoop;      }
     float  meanHitRadialDist() const { return _meanHitRadialDist; }
     float  d0               () const { return _d0;                }
-    float  nstrawhits       () const { return _nStrawHits;        }
+    int    nstrawhits       () const { return _nStrawHits;        }
 
-    //function that evaluates the ratio between the measured ComboHits and the
-    //expected intesections of the helix with the tracker planes. This function
-    //models the tracker as a perfect cylinder
-    float  hitRatio         () const { return _hitRatio;          }
+    //function that evaluates the ratio between the measured straw hits and the
+    //expected intersections of the helix with the tracker plane faces inside the tracker
+    //active radius. If a StrawNoSignal function is given, each face crossing is weighted by
+    //the fraction of live straws (nearest straw in each layer) at the crossing point, so that
+    //dead panels and straws don't lower the ratio. Evaluated on the first call.
+    float  hitRatio         () const;
 
     void   dirOfProp(float& slope, float& slopeErr, float& chi2ndof);
 
   private:
+    float  liveFraction(const Plane& Pln, uint16_t Face, float X, float Y) const;
+
     const HelixSeed* _hel;
-    int        _nMinHitsLoop;
 
     float      _nLoops;
     int        _nStrawHits;
     float      _meanHitRadialDist;
     float      _d0;
-    float      _hitRatio;
-    float      _trackerRIn, _trackerROut, _trackerLength;
+    mutable float _hitRatio = -1.; // < 0: not evaluated yet
+    float      _trackerRIn, _trackerROut;
     const mu2e::Tracker*_tracker;
+    StrawNoSignal _noSignal;
   };
 
 } // namespace mu2e
diff --git a/Mu2eUtilities/src/HelixTool.cc b/Mu2eUtilities/src/HelixTool.cc
index dcfb88070c..e33982dd25 100644
--- a/Mu2eUtilities/src/HelixTool.cc
+++ b/Mu2eUtilities/src/HelixTool.cc
@@ -1,6 +1,8 @@
 #include <cmath>
 #include <stddef.h>
+#include <algorithm>
 #include <array>
+#include <utility>
 #include <vector>
 
 #include "CLHEP/Vector/ThreeVector.h"
@@ -22,9 +24,10 @@
 
 namespace mu2e {
 
-  HelixTool::HelixTool(const HelixSeed *Helix, const mu2e::Tracker*MyTracker) {
+  HelixTool::HelixTool(const HelixSeed *Helix, const mu2e::Tracker*MyTracker, StrawNoSignal NoSignal) {
     _hel           = Helix;
     _tracker       = MyTracker;
+    _noSignal      = std::move(NoSignal);
     _trackerRIn    = _tracker->g4Tracker()->getInnerTrackerEnvelopeParams().innerRadius();
     _trackerROut   = _tracker->g4Tracker()->getInnerTrackerEnvelopeParams().outerRadius();
 
@@ -66,38 +69,97 @@ namespace mu2e {
 
     //here we evaluate once the impact parameter
     _d0 = robustHel->rcent  () - robustHel->radius ();
+  }
+
 
+  float  HelixTool::hitRatio() const {
+    if (_hitRatio >= 0.) return _hitRatio;
 
-    // we now estiamte the ratio of the number of measured hits to the number of the expected ones
-    // we make a few assumptions and appriximations:
+    // ratio of the number of measured hits to the number of expected ones: count the plane
+    // faces the helix crosses inside the tracker active radius, weighted by the fraction of
+    // live straws at the crossing
+    const mu2e::RobustHelix& robustHel = _hel->helix();
     float    expected_faces(1e-6);
     for (size_t planeId=0; planeId<_tracker->nPlanes(); planeId++) {
-      const Plane* pln = &_tracker->getPlane(planeId);
-      int   nPanels = pln->nPanels();
-      if (nPanels == 0 )         continue;
-      std::array<int,2>   idPanels = {0, (nPanels-1)};
+      const Plane& pln = _tracker->getPlane(planeId);
+      if (pln.nPanels() == 0 )         continue;
 
-      for (size_t ipn=0; ipn<idPanels.size(); ++ipn){
-        const Panel* panel = &pln->getPanel(ipn);
-        float    z = (panel->getStraw(0).getMidPoint().z()+panel->getStraw(1).getMidPoint().z())/2.;
-        // if (z < z_first_hit )  continue;
-        // if (z > z_last_hit  )  continue;
+      for (uint16_t face=0; face<StrawId::_nfaces; ++face){
+        // panel number 'face' is on that face
+        const Panel& panel = pln.getPanel(face);
+        float    z = (panel.getStraw(0).getMidPoint().z()+panel.getStraw(1).getMidPoint().z())/2.;
 
         XYZVectorF  pos;
         pos.SetZ(z);
-        robustHel->position(pos);
+        robustHel.position(pos);
 
         //now check that we are in the active area of tracker
         float   hitR = sqrtf(pos.x()*pos.x() + pos.y()*pos.y());
         if (hitR < _trackerRIn )  continue;
         if (hitR > _trackerROut)  continue;
 
-        ++expected_faces;
-      }//end loop over the two panels
+        expected_faces += _noSignal ? liveFraction(pln, face, pos.x(), pos.y()) : 1.;
+      }//end loop over the two faces
     }//end loop over the planes
 
-    //each plane has two layers of panels. We are assuming 2 ComboHits per face
     _hitRatio = _nStrawHits/expected_faces;
+    return _hitRatio;
+  }
+
+
+  float  HelixTool::liveFraction(const Plane& Pln, uint16_t Face, float X, float Y) const {
+    // find the panel of this face covering the crossing point: the one whose outward direction,
+    // perpendicular to its straws in the transverse plane, is closest to the crossing direction
+    const Panel* best(nullptr);
+    float        bestRho(0), nx(0), ny(0);
+    for (int ipn=Face; ipn<int(Pln.nPanels()); ipn+=StrawId::_nfaces) {
+      const Panel& panel = Pln.getPanel(ipn);
+      const auto&  mid   = panel.getStraw(0).getMidPoint();
+      const auto&  dir   = panel.getStraw(0).getDirection();
+      float        along = mid.x()*dir.x() + mid.y()*dir.y();
+      float        px    = mid.x() - along*dir.x();
+      float        py    = mid.y() - along*dir.y();
+      float        norm  = sqrtf(px*px + py*py);
+      if (norm <= 0.)  continue;
+      px /= norm;
+      py /= norm;
+      float        rho   = X*px + Y*py;
+      if (best == nullptr || rho > bestRho) {
+        best    = &panel;
+        bestRho = rho;
+        nx      = px;
+        ny      = py;
+      }
+    }
+    if (best == nullptr) return 1.;
+
+    // straws are ordered radially (straw 0 innermost), alternating between the two layers:
+    // estimate the straw at the crossing and take the nearest straw of each layer around it
+    auto  strawRho = [&](int s) {
+      const auto& mid = best->getStraw(s).getMidPoint();
+      return float(mid.x()*nx + mid.y()*ny);
+    };
+    const int nStraws = best->nStraws();
+    float rho0   = strawRho(0);
+    float rhoMax = strawRho(nStraws-1);
+    int   sEst   = (rhoMax > rho0) ? std::lround((bestRho-rho0)/(rhoMax-rho0)*(nStraws-1)) : 0;
+    sEst = std::clamp(sEst, 0, nStraws-1);
+
+    float nLive(0);
+    for (int layer=0; layer<StrawId::_nlayers; ++layer) {
+      int   nearest(-1);
+      float dmin(0);
+      for (int s=std::max(0, sEst-2); s<=std::min(nStraws-1, sEst+2); ++s) {
+        if (best->getStraw(s).id().layer() != layer) continue;
+        float d = std::fabs(strawRho(s) - bestRho);
+        if (nearest < 0 || d < dmin) {
+          nearest = s;
+          dmin    = d;
+        }
+      }
+      if (nearest >= 0 && !_noSignal(best->getStraw(nearest).id())) ++nLive;
+    }
+    return nLive/StrawId::_nlayers;
   }
 
 
diff --git a/Trigger/CMakeLists.txt b/Trigger/CMakeLists.txt
index 1c8df531b2..28ba03d424 100644
--- a/Trigger/CMakeLists.txt
+++ b/Trigger/CMakeLists.txt
@@ -32,7 +32,9 @@ cet_build_plugin(ReadTriggerInfo art::module
       Offline::GlobalConstantsService
       Offline::MCDataProducts
       Offline::Mu2eUtilities
+      Offline::ProditionsService
       Offline::RecoDataProducts
+      Offline::TrackerConditions
       Offline::TrackerGeom
       KinKal::Trajectory
 )
diff --git a/Trigger/src/ReadTriggerInfo_module.cc b/Trigger/src/ReadTriggerInfo_module.cc
index 6b61db3c49..d69525bcec 100644
--- a/Trigger/src/ReadTriggerInfo_module.cc
+++ b/Trigger/src/ReadTriggerInfo_module.cc
@@ -32,6 +32,8 @@
 #include "Offline/GeometryService/inc/GeomHandle.hh"
 #include "Offline/GeometryService/inc/DetectorSystem.hh"
 #include "Offline/TrackerGeom/inc/Tracker.hh"
+#include "Offline/ProditionsService/inc/ProditionsHandle.hh"
+#include "Offline/TrackerConditions/inc/TrackerStatus.hh"
 
 //Data products
 #include "Offline/RecoDataProducts/inc/CaloCluster.hh"
@@ -260,6 +262,8 @@ namespace mu2e {
     occupancyHist_            _occupancyHist;
 
     const mu2e::Tracker*      _tracker = nullptr;
+    ProditionsHandle<TrackerStatus> _trackerStatus_h;
+    HelixTool::StrawNoSignal  _noSignal;  // straws not counted as expected hits in the helix hit ratio
 
     float  _minPOT, _maxPOT;
     const int _genOccIndex;
@@ -765,6 +769,9 @@ namespace mu2e {
 
   //--------------------------------------------------------------------------------//
   void ReadTriggerInfo::analyze(const art::Event& event) {
+    // same definition of the helix hit ratio as in HelixFilter
+    const TrackerStatus* status = &_trackerStatus_h.get(event.id());
+    _noSignal = [status](StrawId const& sid) { return status->noSignal(sid) || status->suppress(sid); };
 
     if(!_useNGen) _nProcess += 1;
 
@@ -1011,7 +1018,7 @@ namespace mu2e {
   //--------------------------------------------------------------------------------//
   void ReadTriggerInfo::fillHelixTrigInfo(int HelTrigIndex, const HelixSeed* HSeed, helixInfoHist_& Hist) {
     if(!Hist._hHelInfo[HelTrigIndex][0]) throw cet::exception("BADCONFIG") << __func__ << ": Helix histogram index " << HelTrigIndex << " not initialized\n";
-    HelixTool helTool(HSeed, _tracker);
+    HelixTool helTool(HSeed, _tracker, _noSignal);
 
     int        nch       = (int)HSeed->hits().size();
     int        nsh(0);
diff --git a/Trigger/src/SConscript b/Trigger/src/SConscript
index 65b2a77290..f197921806 100644
--- a/Trigger/src/SConscript
+++ b/Trigger/src/SConscript
@@ -23,6 +23,7 @@ helper.make_plugins( [ mainlib,
                        'mu2e_MCDataProducts',
                        'mu2e_RecoDataProducts',
                        'mu2e_TrackerConditions',
+                       'mu2e_DbTables',
                        'mu2e_GeometryService',
                        'mu2e_BFieldGeom',
                        'mu2e_CalorimeterGeom',
diff --git a/TrkFilters/CMakeLists.txt b/TrkFilters/CMakeLists.txt
index d7dbca9c1c..90886f149c 100644
--- a/TrkFilters/CMakeLists.txt
+++ b/TrkFilters/CMakeLists.txt
@@ -6,7 +6,9 @@ cet_build_plugin(HelixFilter art::module
       Offline::GeneralUtilities
       Offline::GeometryService
       Offline::Mu2eUtilities
+      Offline::ProditionsService
       Offline::RecoDataProducts
+      Offline::TrackerConditions
       Offline::TrackerGeom
 )
 
@@ -24,7 +26,9 @@ cet_build_plugin(MultiHelixFilter art::module
       Offline::GeneralUtilities
       Offline::GeometryService
       Offline::Mu2eUtilities
+      Offline::ProditionsService
       Offline::RecoDataProducts
+      Offline::TrackerConditions
       Offline::TrackerGeom
 )
 
diff --git a/TrkFilters/src/HelixFilter_module.cc b/TrkFilters/src/HelixFilter_module.cc
index 4da6f163dd..6606338814 100644
--- a/TrkFilters/src/HelixFilter_module.cc
+++ b/TrkFilters/src/HelixFilter_module.cc
@@ -17,6 +17,8 @@
 #include "Offline/GeometryService/inc/GeomHandle.hh"
 #include "Offline/GeometryService/inc/DetectorSystem.hh"
 #include "Offline/TrackerGeom/inc/Tracker.hh"
+#include "Offline/ProditionsService/inc/ProditionsHandle.hh"
+#include "Offline/TrackerConditions/inc/TrackerStatus.hh"
 // data
 #include "Offline/RecoDataProducts/inc/HelixSeed.hh"
 #include "Offline/DataProducts/inc/Helicity.hh"
@@ -104,6 +106,8 @@ namespace mu2e
       HelixCutsTool():_configured(false){};
 
       void setTrackerGeomHandle(const Tracker* TrackerGeom) { _myTracker = TrackerGeom; }
+      // straws without signal are not counted as expected hits in the hit ratio; nullptr: all straws live
+      void setTrackerStatus(const TrackerStatus* Status) { _trackerStatus = Status; }
 
       int   evalIPAPresc(const float &phi0){
         int val= (_prescalerPar._amplitude - (_prescalerPar._amplitude-1)*sin(_prescalerPar._frequency*phi0 + _prescalerPar._phase));
@@ -116,7 +120,11 @@ namespace mu2e
         //check the helicity
         if (_doHelicityCheck && !(Helix.helix().helicity() == Helicity(_hel)))  return false;
         float mm2MeV = 3./10.;//FIXME!
-        HelixTool helTool(&Helix, _myTracker);
+        HelixTool::StrawNoSignal noSignal;
+        if (_trackerStatus) {
+          noSignal = [status = _trackerStatus](StrawId const& sid) { return status->noSignal(sid) || status->suppress(sid); };
+        }
+        HelixTool helTool(&Helix, _myTracker, noSignal);
         // compute the helix momentum.  Note this is in units of mm!!!
         float hmom       = Helix.helix().momentum()*mm2MeV;
         int   nstrawhits = helTool.nstrawhits();
@@ -126,7 +134,6 @@ namespace mu2e
         float d0         = Helix.helix().rcent() - Helix.helix().radius();
         float lambda     = std::fabs(Helix.helix().lambda());
         float nLoops     = helTool.nLoops();
-        float hRatio     = helTool.hitRatio();
         const float slope          = Helix.recoDir().slope();
         const float slopeErr       = std::fabs(Helix.recoDir().slopeErr());
         const float slopeSignedSig = (slopeErr > 0.f) ? slope/slopeErr : 0.f; //signifance from 0 signed by the slope direction
@@ -135,7 +142,7 @@ namespace mu2e
           std::cout << "[HelixFilter] : status = " << Helix.status() << " nhits = " << nstrawhits << " mom = " << hmom << std::endl
                     << "[HelixFilter] : chi2XY = " << chi2XY << " chi2ZPHI = " << chi2PhiZ << " d0 = " << d0 << " lambda = "<< lambda
                     << " p = " << hmom << " pT = " << hpT
-                    << " nLoops = " << nLoops << " hRatio = "<< hRatio << " slopeSignedSig = " << slopeSignedSig << std::endl;
+                    << " nLoops = " << nLoops << " hRatio = "<< helTool.hitRatio() << " slopeSignedSig = " << slopeSignedSig << std::endl;
         }
         if( Helix.status().hasAllProperties(_goodh)      &&
             (!_hascc || Helix.caloCluster().isNonnull()) &&
@@ -153,7 +160,7 @@ namespace mu2e
             hmom       <= _maxmom        &&
             (!_useSlopeSigMin || slopeSignedSig > _slopeSigMin) &&
             (!_useSlopeSigMax || slopeSignedSig < _slopeSigMax) &&
-            hRatio     >= _minHitRatio ) {
+            (_minHitRatio <= 0. || helTool.hitRatio() >= _minHitRatio) ) { // hit ratio last: it is the most expensive
           //now check if we want to prescake or not
           if (_prescaleUsingD0Phi) {
             float phiAtD0   = Helix.helix().fcent();
@@ -193,6 +200,7 @@ namespace mu2e
       bool          _prescaleUsingD0Phi = false;
       PhiPrescalingParams     _prescalerPar;
       const Tracker*_myTracker     = nullptr;
+      const TrackerStatus* _trackerStatus = nullptr;
     };
 
     struct Config{
@@ -205,6 +213,7 @@ namespace mu2e
       fhicl::Atom<bool>               noFilter             { Name("noFilter"),                Comment("don't apply the filter decision"), false};
       fhicl::Atom<unsigned>           minNHelices          { Name("minNHelices"),             Comment("minimum number of helices passing the cuts"), 1};
       fhicl::Atom<float>              maxDt0               { Name("maxDt0"),                  Comment("maximum difference in time between helices"), -1};
+      fhicl::Atom<bool>               useTrackerStatus     { Name("useTrackerStatus"),        Comment("don't count straws without signal (TrackerStatus) as expected hits in the hit ratio"), true};
     };
 
 
@@ -229,6 +238,8 @@ namespace mu2e
     bool          _noFilter;
     unsigned      _minNHelices;
     float         _maxDt0;
+    bool          _useTrackerStatus;
+    ProditionsHandle<TrackerStatus> _trackerStatus_h;
 
     bool  checkHelixFromHelicity(const HelixSeed&helix);
   };
@@ -243,7 +254,8 @@ namespace mu2e
     _npass(0),
     _noFilter          (config().noFilter()),
     _minNHelices       (config().minNHelices()),
-    _maxDt0            (config().maxDt0())
+    _maxDt0            (config().maxDt0()),
+    _useTrackerStatus  (config().useTrackerStatus())
     {
       produces<TriggerInfo>();
       if(_minNHelices < 2 && _maxDt0 >= 0.) throw cet::exception("BADCONFIG") << "Requested a timing difference cut of " << _maxDt0 << " ns between helices but only " << _minNHelices << " helices required";
@@ -291,6 +303,11 @@ namespace mu2e
     // create output
     std::unique_ptr<TriggerInfo> triginfo(new TriggerInfo);
     ++_nevt;
+    if (_useTrackerStatus) {
+      const TrackerStatus& trackerStatus = _trackerStatus_h.get(evt.id());
+      _posHelCuts.setTrackerStatus(&trackerStatus);
+      _negHelCuts.setTrackerStatus(&trackerStatus);
+    }
     // find the collection
     auto hsH = evt.getValidHandle<HelixSeedCollection>(_hsTag);
     const HelixSeedCollection* hscol = hsH.product();
diff --git a/TrkFilters/src/MultiHelixFilter_module.cc b/TrkFilters/src/MultiHelixFilter_module.cc
index d313cc0cdb..11ba3ab223 100644
--- a/TrkFilters/src/MultiHelixFilter_module.cc
+++ b/TrkFilters/src/MultiHelixFilter_module.cc
@@ -17,6 +17,8 @@
 #include "Offline/GeometryService/inc/GeomHandle.hh"
 #include "Offline/GeometryService/inc/DetectorSystem.hh"
 #include "Offline/TrackerGeom/inc/Tracker.hh"
+#include "Offline/ProditionsService/inc/ProditionsHandle.hh"
+#include "Offline/TrackerConditions/inc/TrackerStatus.hh"
 
 // data
 #include "Offline/RecoDataProducts/inc/CaloCluster.hh"
@@ -104,6 +106,7 @@ namespace mu2e
       }
       void setB(const float bz0) { _bz0 = bz0; }
       void setTracker(const Tracker* tracker) { _tracker = tracker; }
+      void setNoSignal(const HelixTool::StrawNoSignal& noSignal) { _noSignal = noSignal; }
       bool applyCuts(const HelixSeed& helix) {
         if(_debugLevel > 2) std::cout << "Testing helix:\n";
 
@@ -121,7 +124,7 @@ namespace mu2e
         if(_minPt       > 0. && pt  < _minPt)                        return false;
 
         // initialize the helix tool
-        HelixTool hTool(&helix, _tracker);
+        HelixTool hTool(&helix, _tracker, _noSignal);
 
         // geometric selections
         const Helicity hel         = helix.helix().helicity();
@@ -148,16 +151,15 @@ namespace mu2e
         const float chi2XY         = helix.helix().chi2dXY();
         const float chi2PhiZ       = helix.helix().chi2dZPhi();
         const int   nhits          = hTool.nstrawhits();
-        const float hRatio         = hTool.hitRatio();
         const bool  hasCluster     = helix.caloCluster().isNonnull();
         if(_debugLevel > 2) std::cout << " chi2XY = " << chi2XY << " chi2PhiZ = " << chi2PhiZ
-                                      << " nhits = " << nhits << " hitRatio = " << hRatio << " hasCC = " << hasCluster << std::endl;
+                                      << " nhits = " << nhits << " hitRatio = " << hTool.hitRatio() << " hasCC = " << hasCluster << std::endl;
 
         if(chi2XY > _maxChi2XY)                                      return false;
         if(chi2PhiZ > _maxChi2PhiZ)                                  return false;
         if(nhits < _minNStrawHits)                                   return false;
-        if(hRatio < _minHitRatio)                                    return false;
         if(_requireCaloCluster && !hasCluster)                       return false;
+        if(_minHitRatio > 0. && hTool.hitRatio() < _minHitRatio)     return false; // last: the most expensive
 
         if(_debugLevel > 2) std::cout << "  --> passed all cuts\n";
         // passes all cuts
@@ -186,6 +188,7 @@ namespace mu2e
 
       double _bz0                      =  1.   ;
       const Tracker* _tracker          = nullptr;
+      HelixTool::StrawNoSignal _noSignal;          // empty: all straws live
       int    _debugLevel               =  0    ;
     };
 
@@ -202,6 +205,7 @@ namespace mu2e
       fhicl::Atom<float>                             minArcGap            { Name("minArcGap"),           Comment("Minimum distance between helix arcs of same curve direction"), -1.};
       fhicl::Sequence<fhicl::Table<HelixCutsConfig>> HelixCuts            { Name("HelixCuts"),           Comment("Helix cut sets")};
       fhicl::Atom<bool>                              doHistograms         { Name("doHistograms"),        Comment("Make debug histograms"), false};
+      fhicl::Atom<bool>                              useTrackerStatus     { Name("useTrackerStatus"),    Comment("don't count straws without signal (TrackerStatus) as expected hits in the hit ratio"), true};
     };
 
 
@@ -258,9 +262,12 @@ namespace mu2e
     // cut sets
     std::vector<HelixCuts> _helixCuts;
     bool          _doHistograms;
+    bool          _useTrackerStatus;
 
     double         _bz0;
     const Tracker* _tracker;
+    ProditionsHandle<TrackerStatus> _trackerStatus_h;
+    HelixTool::StrawNoSignal        _noSignal; // empty: all straws live
 
     // counters
     unsigned      _nevt, _npass;
@@ -280,6 +287,7 @@ namespace mu2e
     , _maxOverlap        (config().maxOverlap())
     , _minArcGap         (config().minArcGap())
     , _doHistograms      (config().doHistograms())
+    , _useTrackerStatus  (config().useTrackerStatus())
     , _nevt              (0)
     , _npass             (0)
     {
@@ -353,7 +361,7 @@ namespace mu2e
     float helSum = 0.;
     for(auto index : indices) {
       const HelixSeed& helix = hcol->at(index);
-      HelixTool hTool(&helix, _tracker);
+      HelixTool hTool(&helix, _tracker, _noSignal);
       constexpr float mmToMeV = 3./10.;
 
       Hist->p              ->Fill(_bz0 * mmToMeV * helix.helix().momentum());
@@ -516,6 +524,13 @@ namespace mu2e
   bool MultiHelixFilter::filter(art::Event& event){
     ++_nevt;
 
+    // straws without signal are not counted as expected hits in the helix hit ratio
+    if(_useTrackerStatus) {
+      const TrackerStatus* status = &_trackerStatus_h.get(event.id());
+      _noSignal = [status](StrawId const& sid) { return status->noSignal(sid) || status->suppress(sid); };
+      for(auto& cuts : _helixCuts) cuts.setNoSignal(_noSignal);
+    }
+
     // create the output
     std::unique_ptr<TriggerInfo> triginfo(new TriggerInfo);
 
diff --git a/TrkFilters/src/SConscript b/TrkFilters/src/SConscript
index 2cf1887502..0b961dd325 100644
--- a/TrkFilters/src/SConscript
+++ b/TrkFilters/src/SConscript
@@ -16,6 +16,8 @@ mainlib = helper.make_mainlib ( [] )
 # Fixme: split into link lists for each module.
 helper.make_plugins( [ mainlib,
   'mu2e_Mu2eUtilities',
+  'mu2e_TrackerConditions',
+  'mu2e_DbTables',
   'mu2e_SeedService',
   'mu2e_GeometryService',
   'mu2e_BFieldGeom',

From 752460d6ec5873f88d15dee9664bd7d47ea9e1e4 Mon Sep 17 00:00:00 2001
From: michaelmackenzie <michaelmackenzie@u.northwestern.edu>
Date: Wed, 14 Oct 2026 11:45:37 -0500
Subject: [PATCH 2/5] AgnosticHelixFinder: best-of-K candidates, 2pi alias
 test, target-constrained refit

Helix finding options for trackers with fewer hits per track (e.g. the commissioning
tracker, where a CE crosses the active straws once per turn and its hits come in short
groups ~2 pi lambda apart). All default to the previous behavior.

- maxHelixCandidates (1), acceptHelixStrawHits (0): keep the best (most straw hits, then
  the lowest circle + line chi2/dof) of up to K viable candidates per search instead of
  the first; a candidate with >= acceptHelixStrawHits straw hits is saved at once.
- minAliasGap (0 = off): refit the phi-z line with the hits beyond the largest z gap of
  the helix hits shifted by -+2pi; keep the lowest chi2/dof within the momentum range,
  then recover points again.
- targetRefitSigma (0 = off), targetRefitPhiWindow (0.5 rad), targetRefitDphiDzBins (100):
  if the best of the K candidates has |d0| > targetRefitSigma, refit its circle through
  the stopping target (error targetRefitSigma) and find the phi-z line with a dphi/dz scan
  on that circle; the refit replaces the candidate if it has at least as many straw hits.
  It addresses the half-radius / half-lambda solutions that fit the hits of short arcs.

With the Sim_commissioning tracker conditions (and tuned trigger cuts), K = 3,
acceptHelixStrawHits 20, minAliasGap 100 and targetRefitSigma 100 raise the
apr_TrkDe_75m70p_D0200 efficiency for CeM from 52.9% to 57.5% (K + alias) and 61.9%
(with looser phi-z cuts and the target refit), CeP 28.5% -> 39.1%, with no increase of
the NoPrimary D0200 rate; nominal tracker +1.5%. APR finder time 0.036 -> 0.067 ms/event.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 CalPatRec/src/AgnosticHelixFinder_module.cc | 291 +++++++++++++++++++-
 1 file changed, 290 insertions(+), 1 deletion(-)

diff --git a/CalPatRec/src/AgnosticHelixFinder_module.cc b/CalPatRec/src/AgnosticHelixFinder_module.cc
index 38a6c5dbca..f51eaaa653 100644
--- a/CalPatRec/src/AgnosticHelixFinder_module.cc
+++ b/CalPatRec/src/AgnosticHelixFinder_module.cc
@@ -94,6 +94,12 @@ namespace mu2e {
       fhicl::Atom<float>           caloClusterSigma       {Name("caloClusterSigma"     ), Comment("error assigned to calo clust")  };
       fhicl::Atom<int>             minNHelixStrawHits     {Name("minNHelixStrawHits"   ), Comment("straw hit save threshold"    )  };
       fhicl::Atom<int>             minNHelixComboHits     {Name("minNHelixComboHits"   ), Comment("combo hit save threshold"    )  };
+      fhicl::Atom<int>             maxHelixCandidates     {Name("maxHelixCandidates"   ), Comment("save the best (most straw hits) of up to this many candidates per search; 1: save the first"), 1};
+      fhicl::Atom<int>             acceptHelixStrawHits   {Name("acceptHelixStrawHits" ), Comment("save a candidate with at least this many straw hits without searching further; <= 0: never"), 0};
+      fhicl::Atom<float>           minAliasGap            {Name("minAliasGap"          ), Comment("test the 2pi phi-z alias across the largest z gap of a helix's hits if it is at least this (mm); <= 0: off"), 0.f};
+      fhicl::Atom<float>           targetRefitSigma       {Name("targetRefitSigma"     ), Comment("if the best of the maxHelixCandidates (> 1) candidates has |d0| > this (mm), refit its circle through the stopping target (error = this) with a dphi/dz scan; the refit replaces it if it has at least as many straw hits; <= 0: off"), 0.f};
+      fhicl::Atom<float>           targetRefitPhiWindow   {Name("targetRefitPhiWindow" ), Comment("phi window (rad) for hits to count on a line in the target refit's dphi/dz scan"), 0.5f};
+      fhicl::Atom<int>             targetRefitDphiDzBins  {Name("targetRefitDphiDzBins"), Comment("dphi/dz bins per helicity in the target refit's dphi/dz scan"), 100};
       fhicl::Atom<float>           minHelixPerpMomentum   {Name("minHelixPerpMomentum" ), Comment("min pt of helix"             )  };
       fhicl::Atom<float>           maxHelixPerpMomentum   {Name("maxHelixPerpMomentum" ), Comment("max pt of helix"             )  };
       fhicl::Atom<float>           minHelixMomentum       {Name("minHelixMomentum"     ), Comment("min momentum of helix"       )  };
@@ -174,6 +180,12 @@ namespace mu2e {
     float    _caloClusterSigmaSq;
     int      _minNHelixStrawHits;
     int      _minNHelixComboHits;
+    int      _maxHelixCandidates;
+    int      _acceptHelixStrawHits;
+    float    _minAliasGap;
+    float    _targetRefitSigma;
+    float    _targetRefitPhiWindow;
+    int      _targetRefitDphiDzBins;
     float    _minHelixPerpMomentum;
     float    _maxHelixPerpMomentum;
     float    _minHelixMomentum;
@@ -211,6 +223,17 @@ namespace mu2e {
     bool                          _intenseEvent;
     bool                          _intenseCluster;
 
+    // best helix candidate of the current search (maxHelixCandidates > 1)
+    struct helixCandidate {
+      std::vector<cHit> tcHits;
+      ::LsqSums4 circleFitter;
+      ::LsqSums2 lineFitter;
+      int nStrawHits = 0, nComboHits = 0, nStrawHitsInTC = 0, nComboHitsInTC = 0;
+      float chi2 = 0.f;
+      float xc = 0.f, yc = 0.f, r = 0.f; // circle
+    };
+    std::vector<helixCandidate>   _candidates;
+
     //-----------------------------------------------------------------------------
     // stuff for tool
     //-----------------------------------------------------------------------------
@@ -264,8 +287,14 @@ namespace mu2e {
     bool         refinePhiLine             (size_t lineIndex);
     void         initFinalSeed             (LoopCondition& outcome);
     bool         recoverPoints             ();
+    bool         resolveGapAlias           ();
+    bool         scanDphiDz                (float window);
+    bool         targetRefit               (const helixCandidate& best, helixCandidate& refit);
     void         checkHelixViability       (LoopCondition& outcome);
     void         saveHelix                 (size_t tc, HelixSeedCollection& HSColl);
+    bool         saveAndContinue           (size_t tc, HelixSeedCollection& HSColl, int nStrawHitsInHelix, int nComboHitsInHelix,
+                                            int nStrawHitsInTimeCluster, int nComboHitsInTimeCluster);
+    bool         saveCandidates            (size_t tc, HelixSeedCollection& HSColl);
     bool         validHelixDirection       (TrkFitDirection::FitDirection direction);
     bool         goodTimeCluster           (const TimeCluster& tc);
 
@@ -317,6 +346,12 @@ namespace mu2e {
     _caloClusterSigma              (config().caloClusterSigma()                      ),
     _minNHelixStrawHits            (config().minNHelixStrawHits()                    ),
     _minNHelixComboHits            (config().minNHelixComboHits()                    ),
+    _maxHelixCandidates            (config().maxHelixCandidates()                    ),
+    _acceptHelixStrawHits          (config().acceptHelixStrawHits()                  ),
+    _minAliasGap                   (config().minAliasGap()                           ),
+    _targetRefitSigma              (config().targetRefitSigma()                      ),
+    _targetRefitPhiWindow          (config().targetRefitPhiWindow()                  ),
+    _targetRefitDphiDzBins         (config().targetRefitDphiDzBins()                 ),
     _minHelixPerpMomentum          (config().minHelixPerpMomentum()                  ),
     _maxHelixPerpMomentum          (config().maxHelixPerpMomentum()                  ),
     _minHelixMomentum              (config().minHelixMomentum()                      ),
@@ -824,6 +859,7 @@ namespace mu2e {
     // now we loop over triplets
     LoopCondition loopCondition;
     triplet tripletInfo;
+    _candidates.clear(); // viable helix candidates of this search (maxHelixCandidates > 1)
     for (size_t i = 0; i < _tcHits.size() - 2; i++) {
       if(_doTiming > 1) _watch->Increment("triplet-i");
       if(_doTiming > 2) _watch->StopTime ("triplet-j");
@@ -868,6 +904,7 @@ namespace mu2e {
           initFinalSeed(loopCondition);
           if (loopCondition == CONTINUE) { continue; }
           while (recoverPoints()); // while points are recovered, keep trying to recover
+          if (_minAliasGap > 0.f && resolveGapAlias()) { while (recoverPoints()); }
           checkHelixViability(loopCondition);
           if (loopCondition == CONTINUE) { continue; }
           // before saving helix we make sure it has enough hits
@@ -887,7 +924,28 @@ namespace mu2e {
             nStrawHitsInHelix += nStrawHitsInHit;
             ++nComboHitsInHelix;
           }
-          if (nStrawHitsInHelix >= _minNHelixStrawHits && nComboHitsInHelix >= _minNHelixComboHits) {
+          if (nStrawHitsInHelix >= _minNHelixStrawHits && nComboHitsInHelix >= _minNHelixComboHits &&
+              _maxHelixCandidates > 1 && !(_acceptHelixStrawHits > 0 && nStrawHitsInHelix >= _acceptHelixStrawHits)) {
+            // keep the best candidate (most straw hits, then the best fit) and continue the search
+            helixCandidate c;
+            c.tcHits = _tcHits;
+            c.circleFitter = _circleFitter;
+            c.lineFitter = _lineFitter;
+            c.nStrawHits = nStrawHitsInHelix;
+            c.nComboHits = nComboHitsInHelix;
+            c.nStrawHitsInTC = nStrawHitsInTimeCluster;
+            c.nComboHitsInTC = nComboHitsInTimeCluster;
+            c.chi2 = _circleFitter.chi2DofCircle() + _lineFitter.chi2Dof();
+            c.xc = _circleFitter.x0(); c.yc = _circleFitter.y0(); c.r = _circleFitter.radius();
+            _candidates.push_back(std::move(c));
+            if ((int)_candidates.size() < _maxHelixCandidates) continue;
+            // enough candidates: save the best one
+            if(_doTiming > 0) _watch->StopTime(__func__);
+            if(_doTiming > 1) _watch->StopTime("triplet-i");
+            if(_doTiming > 2) _watch->StopTime("triplet-j");
+            if(_doTiming > 3) _watch->StopTime("triplet-k");
+            return saveCandidates(tc, HSColl);
+          } else if (nStrawHitsInHelix >= _minNHelixStrawHits && nComboHitsInHelix >= _minNHelixComboHits) {
             saveHelix(tc, HSColl);
             if (_diagLevel > 0) { _diagInfo.nHelices++; }
             // we only want to search for another helix if we have enough remaining hits after saving helix
@@ -915,10 +973,102 @@ namespace mu2e {
     if(_doTiming > 2) _watch->StopTime("triplet-j");
     if(_doTiming > 3) _watch->StopTime("triplet-k");
 
+    // search finished with fewer than maxHelixCandidates candidates: save the best one
+    if (!_candidates.empty()) return saveCandidates(tc, HSColl);
+
     // No helix found, no need to continue searching in this time cluster
     return false;
   }
 
+  //-----------------------------------------------------------------------------
+  // target-constrained refit of a candidate with |d0| > targetRefitSigma: the circle is refit
+  // with its hits plus the stopping target (0,0) with error targetRefitSigma, then the phi-z line
+  // is found by a dphi/dz scan on that circle, the hits recovered and the helix checked as usual.
+  // Short arcs leave the radius and the 2pi phi-z count ambiguous (half radius / half lambda
+  // solutions fit the hits too); this tries the solution through the target.
+  //-----------------------------------------------------------------------------
+  bool AgnosticHelixFinder::targetRefit(const helixCandidate& best, helixCandidate& refit) {
+    const float d0 = std::sqrt(best.xc * best.xc + best.yc * best.yc) - best.r;
+    if (std::fabs(d0) <= _targetRefitSigma) return false;
+    _tcHits = best.tcHits;
+    const float wT = 1.f / (_targetRefitSigma * _targetRefitSigma);
+    // circle: the candidate's hits plus the target, errors updated once
+    for (int iter = 0; iter < 2; ++iter) {
+      const float xC = iter == 0 ? best.xc : _circleFitter.x0();
+      const float yC = iter == 0 ? best.yc : _circleFitter.y0();
+      _circleFitter.clear();
+      _circleFitter.addPoint(0.f, 0.f, wT);
+      for (size_t i = 0; i < _tcHits.size(); i++) {
+        const auto& h = _tcHits[i];
+        if (!h.used || h.inHelix || h.hitIndice == HitType::STOPPINGTARGET || h.hitIndice == HitType::CALOCLUSTER) continue;
+        computeCircleError2(i, xC, yC);
+        _circleFitter.addPoint(getPos(i).x(), getPos(i).y(), 1.f / _tcHits[i].circleError2);
+      }
+    }
+    const float pt = computeHelixPerpMomentum(_circleFitter.radius());
+    if (pt < _minHelixPerpMomentum || pt > _maxHelixPerpMomentum) return false;
+    if (!scanDphiDz(_targetRefitPhiWindow)) return false;
+    while (recoverPoints());
+    LoopCondition outcome;
+    checkHelixViability(outcome);
+    if (outcome == CONTINUE) return false;
+    int nsh = 0, nch = 0;
+    for (const auto& h : _tcHits) {
+      if (!h.used || h.inHelix || h.hitIndice == HitType::STOPPINGTARGET || h.hitIndice == HitType::CALOCLUSTER) continue;
+      nsh += h.hit->nStrawHits();
+      ++nch;
+    }
+    if (nsh < _minNHelixStrawHits || nch < _minNHelixComboHits) return false;
+    refit.tcHits = _tcHits;
+    refit.circleFitter = _circleFitter;
+    refit.lineFitter = _lineFitter;
+    refit.nStrawHits = nsh;
+    refit.nComboHits = nch;
+    refit.nStrawHitsInTC = best.nStrawHitsInTC;
+    refit.nComboHitsInTC = best.nComboHitsInTC;
+    refit.chi2 = _circleFitter.chi2DofCircle() + _lineFitter.chi2Dof();
+    refit.xc = _circleFitter.x0(); refit.yc = _circleFitter.y0(); refit.r = _circleFitter.radius();
+    return true;
+  }
+
+  //-----------------------------------------------------------------------------
+  // save the best candidate (most straw hits, then the best fit), or its target-constrained
+  // refit (targetRefitSigma > 0) if that has at least as many straw hits. Returns whether to
+  // search for another helix.
+  //-----------------------------------------------------------------------------
+  bool AgnosticHelixFinder::saveCandidates(size_t tc, HelixSeedCollection& HSColl) {
+    size_t ib = 0;
+    for (size_t i = 1; i < _candidates.size(); ++i) {
+      const auto& c = _candidates[i];
+      const auto& b = _candidates[ib];
+      if (c.nStrawHits > b.nStrawHits || (c.nStrawHits == b.nStrawHits && c.chi2 < b.chi2)) ib = i;
+    }
+    if (_targetRefitSigma > 0.f) {
+      helixCandidate refit;
+      if (targetRefit(_candidates[ib], refit) && refit.nStrawHits >= _candidates[ib].nStrawHits)
+        _candidates[ib] = std::move(refit);
+    }
+    const auto& best = _candidates[ib];
+    _tcHits = best.tcHits;
+    _circleFitter = best.circleFitter;
+    _lineFitter = best.lineFitter;
+    return saveAndContinue(tc, HSColl, best.nStrawHits, best.nComboHits, best.nStrawHitsInTC, best.nComboHitsInTC);
+  }
+
+  //-----------------------------------------------------------------------------
+  // save the current helix (fitters and used hits) and decide whether to search for another one
+  //-----------------------------------------------------------------------------
+  bool AgnosticHelixFinder::saveAndContinue(size_t tc, HelixSeedCollection& HSColl, int nStrawHitsInHelix, int nComboHitsInHelix,
+                                            int nStrawHitsInTimeCluster, int nComboHitsInTimeCluster) {
+    saveHelix(tc, HSColl);
+    if (_diagLevel > 0) { _diagInfo.nHelices++; }
+    const int remainingStrawHits = nStrawHitsInTimeCluster - nStrawHitsInHelix;
+    const int remainingComboHits = nComboHitsInTimeCluster - nComboHitsInHelix;
+    const bool findAnotherHelix = remainingStrawHits >= _minNHelixStrawHits && remainingComboHits >= _minNHelixComboHits;
+    if(findAnotherHelix && _doAverageFlag) resetFlags(); // reset flags if needed
+    return findAnotherHelix;
+  }
+
   //-----------------------------------------------------------------------------
   // check flags to see if point is good for triplet-ing with
   //-----------------------------------------------------------------------------
@@ -1567,6 +1717,145 @@ namespace mu2e {
     return recoveries;
   }
 
+  //-----------------------------------------------------------------------------
+  // test the 2pi alias of the phi-z line across the largest z gap between the used hits: the
+  // phi of the hits on one side of the gap is shifted by +-2pi and the line refit. The
+  // alternative with the lowest chi2/dof, if lower than the current line's and within the
+  // momentum range, replaces it. Returns true if the line changed.
+  //-----------------------------------------------------------------------------
+  bool AgnosticHelixFinder::resolveGapAlias() {
+    if(_doTiming > 4) _watch->SetTime(__func__);
+    // used tracker hits, in decreasing z (the _tcHits order)
+    std::vector<size_t> idx;
+    for (size_t i = 0; i < _tcHits.size(); i++) {
+      const auto& h = _tcHits[i];
+      if (!h.used || h.inHelix || h.hitIndice == HitType::STOPPINGTARGET || h.hitIndice == HitType::CALOCLUSTER) continue;
+      idx.push_back(i);
+    }
+    if (idx.size() < 4) { if(_doTiming > 4) _watch->StopTime(__func__); return false; }
+    size_t igap = 0; // hits idx[0..igap] are on the high-z side
+    float gap = 0.f;
+    for (size_t k = 0; k + 1 < idx.size(); ++k) {
+      const float dz = getPos(idx[k]).z() - getPos(idx[k+1]).z();
+      if (dz > gap) { gap = dz; igap = k; }
+    }
+    if (gap < _minAliasGap) { if(_doTiming > 4) _watch->StopTime(__func__); return false; }
+
+    const float radius = _circleFitter.radius();
+    float bestChi2 = _lineFitter.chi2Dof();
+    int   bestShift = 0;
+    ::LsqSums2 bestFitter;
+    for (int shift : {-1, 1}) {
+      ::LsqSums2 fitter;
+      for (size_t k = 0; k < idx.size(); ++k) {
+        const auto& h = _tcHits[idx[k]];
+        const int n = h.helixPhiCorrection + (k > igap ? shift : 0);
+        fitter.addPoint(getPos(idx[k]).z(), h.helixPhi + _twopi * n, 1.f / h.helixPhiError2);
+      }
+      const float mom2 = computeHelixMomentum2(radius, fitter.dydx());
+      if (mom2 < _minHelixMomentumSq || mom2 > _maxHelixMomentumSq) continue;
+      if (fitter.chi2Dof() < bestChi2) {
+        bestChi2 = fitter.chi2Dof();
+        bestShift = shift;
+        bestFitter = fitter;
+      }
+    }
+    if (bestShift != 0) {
+      for (size_t k = igap + 1; k < idx.size(); ++k) _tcHits[idx[k]].helixPhiCorrection += bestShift;
+      _lineFitter = bestFitter;
+    }
+    if(_doTiming > 4) _watch->StopTime(__func__);
+    return bestShift != 0;
+  }
+
+  //-----------------------------------------------------------------------------
+  // scan the phi-z slope on the current circle (target refit): for each dphi/dz allowed by the
+  // momentum range (both helicities, targetRefitDphiDzBins bins each), count the time-cluster
+  // hits on the circle (within maxCircleRecoverSigma) whose phi lies on a common line, modulo
+  // 2pi, within window. The hits on the best line replace the current ones and the line and
+  // circle are refit. Returns false if no line with >= 3 hits is found.
+  //-----------------------------------------------------------------------------
+  bool AgnosticHelixFinder::scanDphiDz(float window) {
+    if(_doTiming > 4) _watch->SetTime(__func__);
+    const float xC = _circleFitter.x0(), yC = _circleFitter.y0(), rC = _circleFitter.radius();
+    const float maxRes2 = _maxCircleRecoverSigma * _maxCircleRecoverSigma;
+    std::vector<size_t> cand;
+    for (size_t i = 0; i < _tcHits.size(); i++) {
+      auto& h = _tcHits[i];
+      if (h.inHelix || h.hitIndice == HitType::STOPPINGTARGET || h.hitIndice == HitType::CALOCLUSTER) continue;
+      computeCircleError2(i, xC, yC);
+      if (computeCircleResidual2(i, xC, yC, rC) > maxRes2) continue;
+      cand.push_back(i);
+    }
+    auto done = [&](bool res) { if(_doTiming > 4) _watch->StopTime(__func__); return res; };
+    if (cand.size() < 3) return done(false);
+    std::vector<float> z(cand.size()), phi(cand.size());
+    for (size_t k = 0; k < cand.size(); ++k) {
+      computeHelixPhi(cand[k], xC, yC);
+      z[k] = getPos(cand[k]).z();
+      phi[k] = _tcHits[cand[k]].helixPhi;
+    }
+    // slope range from the momentum limits: p^2 = (R^2 + lambda^2) (cB)^2, |dphi/dz| <= maxDphiDz
+    const float lam2Max = _maxHelixMomentumSq / _bz0ConvSq - rC * rC;
+    const float lam2Min = _minHelixMomentumSq / _bz0ConvSq - rC * rC;
+    const float sLow  = (lam2Max > 0.f) ? 1.f / std::sqrt(lam2Max) : 0.f;
+    float sHigh = _maxDphiDz;
+    if (lam2Min > 0.f) sHigh = std::min(sHigh, 1.f / std::sqrt(lam2Min));
+    if (sHigh <= sLow) return done(false);
+
+    int bestN = 0;
+    float bestS = 0.f, bestPhi0 = 0.f;
+    std::vector<float> ph(cand.size()), ph2(2 * cand.size());
+    const int nb = std::max(1, _targetRefitDphiDzBins);
+    for (int sign : {-1, 1}) {
+      for (int b = 0; b <= nb; ++b) {
+        const float slope = sign * (sLow + (sHigh - sLow) * b / nb);
+        for (size_t k = 0; k < cand.size(); ++k) {
+          float r = std::fmod(phi[k] - slope * z[k], float(_twopi));
+          if (r < 0.f) r += _twopi;
+          ph[k] = r;
+        }
+        std::sort(ph.begin(), ph.end());
+        for (size_t k = 0; k < ph.size(); ++k) { ph2[k] = ph[k]; ph2[k + ph.size()] = ph[k] + _twopi; }
+        // largest number of phases within a window (circular)
+        size_t j = 0;
+        for (size_t k = 0; k < ph.size(); ++k) {
+          if (j < k) j = k;
+          while (j + 1 < ph2.size() && ph2[j + 1] - ph2[k] <= window) ++j;
+          const int n = j - k + 1;
+          if (n > bestN) { bestN = n; bestS = slope; bestPhi0 = 0.5f * (ph2[k] + ph2[j]); }
+        }
+      }
+    }
+    if (bestN < 3) return done(false);
+
+    // adopt the best line: hits within the window are used, with their 2pi corrections
+    for (size_t i = 0; i < _tcHits.size(); i++) {
+      auto& h = _tcHits[i];
+      if (h.inHelix || h.hitIndice == HitType::STOPPINGTARGET || h.hitIndice == HitType::CALOCLUSTER) continue;
+      h.used = false;
+    }
+    _lineFitter.clear();
+    _circleFitter.clear();
+    for (size_t k = 0; k < cand.size(); ++k) {
+      auto& h = _tcHits[cand[k]];
+      const float pred = bestS * z[k] + bestPhi0;
+      const int n = std::round((pred - phi[k]) / _twopi);
+      if (std::fabs(phi[k] + _twopi * n - pred) > 0.5f * window + 1.e-4f) continue;
+      h.used = true;
+      h.helixPhiCorrection = n;
+      _lineFitter.addPoint(z[k], phi[k] + _twopi * n, 1.f / h.helixPhiError2);
+      _circleFitter.addPoint(getPos(cand[k]).x(), getPos(cand[k]).y(), 1.f / h.circleError2);
+    }
+    // keep the calo cluster in the circle fit, as in the other fits
+    for (size_t i = 0; i < _tcHits.size(); i++) {
+      const auto& h = _tcHits[i];
+      if (h.hitIndice == HitType::CALOCLUSTER && h.used && !h.inHelix)
+        _circleFitter.addPoint(getPos(i).x(), getPos(i).y(), 1.f / h.circleError2);
+    }
+    return done(true);
+  }
+
   //-----------------------------------------------------------------------------
   // check if helix is savable based on fcl defined cuts
   //-----------------------------------------------------------------------------

From 3aa9fd8776f40e409366c1aa83c6c6e752844013 Mon Sep 17 00:00:00 2001
From: michaelmackenzie <michaelmackenzie@u.northwestern.edu>
Date: Wed, 14 Oct 2026 11:45:37 -0500
Subject: [PATCH 3/5] Analyses: AprHelixTruthAna, helix seeds vs MC truth of
 the signal particle

Per-event TTree comparing a helix seed collection (e.g. the trigger APR helices) with the
MC truth of the main signal particle (the highest-momentum e+- or mu+- with >= minHits
straw hits): its momentum, direction, circle and lambda at its first straw gas step, its
combo-hit z span and gaps, and for each helix its parameters, straw hits and signal straw
hits, circle chi2 with the reco and the true circle, plus the signal hits themselves.
Optionally records the decisions of trigger paths of the current process (triggerPaths),
so trigger efficiencies can be split by the signal truth (particle type, direction).
Used to study the APR helix finding with the commissioning tracker conditions.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 Analyses/CMakeLists.txt                 |  13 +
 Analyses/src/AprHelixTruthAna_module.cc | 330 ++++++++++++++++++++++++
 2 files changed, 343 insertions(+)
 create mode 100644 Analyses/src/AprHelixTruthAna_module.cc

diff --git a/Analyses/CMakeLists.txt b/Analyses/CMakeLists.txt
index 1520f73fb9..4aedb8b320 100644
--- a/Analyses/CMakeLists.txt
+++ b/Analyses/CMakeLists.txt
@@ -1,3 +1,16 @@
+cet_build_plugin(AprHelixTruthAna art::module
+    REG_SOURCE src/AprHelixTruthAna_module.cc
+    LIBRARIES REG
+      art_root_io::TFileService_service
+
+      Offline::BFieldGeom
+      Offline::GeometryService
+      Offline::MCDataProducts
+      Offline::Mu2eUtilities
+      Offline::RecoDataProducts
+      ROOT::Tree
+)
+
 cet_build_plugin(BFieldPlotter art::module
     REG_SOURCE src/BFieldPlotter_module.cc
     LIBRARIES REG
diff --git a/Analyses/src/AprHelixTruthAna_module.cc b/Analyses/src/AprHelixTruthAna_module.cc
new file mode 100644
index 0000000000..b07b36af77
--- /dev/null
+++ b/Analyses/src/AprHelixTruthAna_module.cc
@@ -0,0 +1,330 @@
+///////////////////////////////////////////////////////////////////////////////
+// AprHelixTruthAna: per-event TTree comparing helix seeds (e.g. the trigger APR
+// helices) with the MC truth of the main signal particle, to study helix finding
+// efficiency and wrong-lambda / wrong-circle helices. Optionally records the
+// decisions of trigger paths of the current process (triggerPaths), so trigger
+// efficiencies can be split by the truth of the signal particle.
+//
+// The signal particle is the highest-momentum e+- or mu+- with >= minHits straw hits
+// (via the StrawDigiMCs). Its truth at the earliest StrawGasStep gives p, pT, pz, the
+// circle center and radius, and lambda = pz/(c B). Its combo hits in the
+// helix finder's input collection give the z span and the largest z gap.
+// For each helix: lambda, radius, center, N(straw hits), N(signal straw hits),
+// z span and largest z gap of its hits.
+///////////////////////////////////////////////////////////////////////////////
+
+#include "art/Framework/Core/EDAnalyzer.h"
+#include "art/Framework/Principal/Event.h"
+#include "art/Framework/Principal/Handle.h"
+#include "art/Framework/Principal/Run.h"
+#include "art/Framework/Services/Registry/ServiceHandle.h"
+#include "art_root_io/TFileService.h"
+#include "fhiclcpp/types/Atom.h"
+
+#include "Offline/BFieldGeom/inc/BFieldManager.hh"
+#include "Offline/GeometryService/inc/DetectorSystem.hh"
+#include "Offline/GeometryService/inc/GeomHandle.hh"
+#include "Offline/MCDataProducts/inc/StrawDigiMC.hh"
+#include "Offline/RecoDataProducts/inc/ComboHit.hh"
+#include "Offline/RecoDataProducts/inc/HelixSeed.hh"
+#include "Offline/Mu2eUtilities/inc/TriggerResultsNavigator.hh"
+#include "canvas/Persistency/Common/TriggerResults.h"
+#include "fhiclcpp/types/Sequence.h"
+
+#include "CLHEP/Units/PhysicalConstants.h"
+#include "TTree.h"
+
+#include <algorithm>
+#include <cmath>
+#include <map>
+#include <vector>
+
+namespace mu2e {
+
+  class AprHelixTruthAna : public art::EDAnalyzer {
+  public:
+    struct Config {
+      using Name = fhicl::Name;
+      using Comment = fhicl::Comment;
+      fhicl::Atom<art::InputTag> hsTag {Name("helixSeedTag"), Comment("helix seed collection")};
+      fhicl::Atom<art::InputTag> chTag {Name("comboHitTag"), Comment("helix finder input combo hits")};
+      fhicl::Atom<art::InputTag> sdmcTag {Name("strawDigiMCTag"), Comment("StrawDigiMC collection")};
+      fhicl::Atom<int> minHits {Name("minHits"), Comment("min straw hits of the signal particle"), 5};
+      fhicl::Sequence<std::string> paths {Name("triggerPaths"), Comment("trigger paths of this process to record (pass = 1)"), std::vector<std::string>{}};
+    };
+    using Parameters = art::EDAnalyzer::Table<Config>;
+    explicit AprHelixTruthAna(const Parameters& config);
+    void beginJob() override;
+    void beginRun(const art::Run&) override;
+    void analyze(const art::Event& event) override;
+
+  private:
+    // z span, largest gap and N(distinct z planes) of a set of z values
+    static void zInfo(std::vector<float> z, float& zmin, float& zmax, float& gap, int& nz);
+
+    art::InputTag _hsTag, _chTag, _sdmcTag;
+    int _minHits;
+    std::vector<std::string> _paths;
+    std::vector<int> _pass; // per entry of _paths
+    float _bzConv = 0.; // mm per MeV/c
+    TTree* _tree = nullptr;
+
+    // event
+    int _run, _subrun, _event;
+    bool _hasHelices; // the helix collection exists (its producer ran)
+    // truth of the signal particle
+    int _pdg, _nsh; // pdg, N(straw hits)
+    float _p, _pt, _pz, _rtrue, _lambdatrue, _xctrue, _yctrue, _d0true;
+    int _nchsig, _nzsig; // signal combo hits in the input collection, distinct z
+    float _zminsig, _zmaxsig, _gapsig;
+    // helices
+    std::vector<float> _hlambda, _hradius, _hxc, _hyc, _hd0, _hzmin, _hzmax, _hgap, _hchi2xy, _hchi2zphi;
+    std::vector<int> _hnsh, _hnshsig, _hnz;
+    // circle chi2/ndof of the helix's hits for the reco and the true circle (finder error model),
+    // and N(signal combo hits in the input collection) within 3 sigma of each circle
+    std::vector<float> _hchi2creco, _hchi2ctrue;
+    std::vector<int> _hnsig3reco, _hnsig3true;
+    int _nsig3true;
+    // signal combo hits: position, N(straw hits), bit mask of the helices containing it
+    std::vector<float> _sx, _sy, _sz;
+    std::vector<int> _snsh, _shel;
+    static float circleChi2(const ComboHit& hit, float xc, float yc, float r);
+  };
+
+  AprHelixTruthAna::AprHelixTruthAna(const Parameters& config) :
+    art::EDAnalyzer{config}, _hsTag(config().hsTag()), _chTag(config().chTag()), _sdmcTag(config().sdmcTag()), _minHits(config().minHits()), _paths(config().paths()), _pass(_paths.size(), 0) {}
+
+  void AprHelixTruthAna::beginJob() {
+    art::ServiceHandle<art::TFileService> tfs;
+    _tree = tfs->make<TTree>("helices", "helix seeds vs signal truth");
+    _tree->Branch("run", &_run);
+    _tree->Branch("subrun", &_subrun);
+    _tree->Branch("event", &_event);
+    _tree->Branch("hasHelices", &_hasHelices);
+    _tree->Branch("pdg", &_pdg);
+    _tree->Branch("nsh", &_nsh);
+    _tree->Branch("p", &_p);
+    _tree->Branch("pt", &_pt);
+    _tree->Branch("pz", &_pz);
+    _tree->Branch("rtrue", &_rtrue);
+    _tree->Branch("lambdatrue", &_lambdatrue);
+    _tree->Branch("xctrue", &_xctrue);
+    _tree->Branch("yctrue", &_yctrue);
+    _tree->Branch("d0true", &_d0true);
+    _tree->Branch("nchsig", &_nchsig);
+    _tree->Branch("nzsig", &_nzsig);
+    _tree->Branch("zminsig", &_zminsig);
+    _tree->Branch("zmaxsig", &_zmaxsig);
+    _tree->Branch("gapsig", &_gapsig);
+    _tree->Branch("hlambda", &_hlambda);
+    _tree->Branch("hradius", &_hradius);
+    _tree->Branch("hxc", &_hxc);
+    _tree->Branch("hyc", &_hyc);
+    _tree->Branch("hd0", &_hd0);
+    _tree->Branch("hzmin", &_hzmin);
+    _tree->Branch("hzmax", &_hzmax);
+    _tree->Branch("hgap", &_hgap);
+    _tree->Branch("hchi2xy", &_hchi2xy);
+    _tree->Branch("hchi2zphi", &_hchi2zphi);
+    _tree->Branch("hnsh", &_hnsh);
+    _tree->Branch("hnshsig", &_hnshsig);
+    _tree->Branch("hnz", &_hnz);
+    _tree->Branch("hchi2creco", &_hchi2creco);
+    _tree->Branch("hchi2ctrue", &_hchi2ctrue);
+    _tree->Branch("hnsig3reco", &_hnsig3reco);
+    _tree->Branch("hnsig3true", &_hnsig3true);
+    _tree->Branch("nsig3true", &_nsig3true);
+    _tree->Branch("sx", &_sx);
+    _tree->Branch("sy", &_sy);
+    _tree->Branch("sz", &_sz);
+    _tree->Branch("snsh", &_snsh);
+    _tree->Branch("shel", &_shel);
+    for (size_t i = 0; i < _paths.size(); ++i) _tree->Branch(("trig_" + _paths[i]).c_str(), &_pass[i]);
+  }
+
+  void AprHelixTruthAna::beginRun(const art::Run&) {
+    GeomHandle<BFieldManager> bfmgr;
+    GeomHandle<DetectorSystem> det;
+    const double bz = bfmgr->getBField(det->toMu2e(CLHEP::Hep3Vector(0., 0., 0.))).z();
+    _bzConv = bz * CLHEP::c_light / 1000.; // MeV/c per mm
+  }
+
+  // squared circle residual in units of the AgnosticHelixFinder circle error
+  float AprHelixTruthAna::circleChi2(const ComboHit& hit, float xc, float yc, float r) {
+    const float dx = hit.pos().x() - xc, dy = hit.pos().y() - yc;
+    const float d2 = dx*dx + dy*dy;
+    const float dxn = dx * hit.vDir().x() + dy * hit.vDir().y();
+    const float costh2 = dxn * dxn / d2;
+    const float err2 = hit.wireVar() * (1.f - costh2) + hit.transVar() * costh2;
+    const float res = r - std::sqrt(d2);
+    return res * res / err2;
+  }
+
+  void AprHelixTruthAna::zInfo(std::vector<float> z, float& zmin, float& zmax, float& gap, int& nz) {
+    zmin = zmax = gap = 0.;
+    nz = 0;
+    if (z.empty()) return;
+    std::sort(z.begin(), z.end());
+    zmin = z.front();
+    zmax = z.back();
+    nz = 1;
+    for (size_t i = 1; i < z.size(); ++i) {
+      const float dz = z[i] - z[i-1];
+      gap = std::max(gap, dz);
+      if (dz > 10.) ++nz; // a new plane face
+    }
+  }
+
+  void AprHelixTruthAna::analyze(const art::Event& event) {
+    _run = event.run(); _subrun = event.subRun(); _event = event.event();
+    if (!_paths.empty()) {
+      art::Handle<art::TriggerResults> trH;
+      event.getByLabel(art::InputTag("TriggerResults", "", moduleDescription().processName()), trH);
+      for (size_t i = 0; i < _paths.size(); ++i) {
+        _pass[i] = 0;
+        if (trH.isValid()) {
+          TriggerResultsNavigator nav(trH.product());
+          if (nav.validPath(_paths[i])) _pass[i] = nav.accepted(_paths[i]);
+        }
+      }
+    }
+    const auto& sdmcs = *event.getValidHandle<StrawDigiMCCollection>(_sdmcTag);
+    art::Handle<ComboHitCollection> chH;
+    event.getByLabel(_chTag, chH);
+    static const ComboHitCollection noHits;
+    const auto& chcol = chH.isValid() ? *chH : noHits;
+    // the helix finder may not have run (an earlier filter in its path failed)
+    art::Handle<HelixSeedCollection> hsH;
+    event.getByLabel(_hsTag, hsH);
+    _hasHelices = hsH.isValid();
+    static const HelixSeedCollection empty;
+    const auto& hscol = _hasHelices ? *hsH : empty;
+
+    // signal particle: the highest-momentum e+- or mu+- (at its earliest straw hit) with >= minHits straw hits
+    // (the conversion electron in CE samples; pileup loopers have more hits but low momentum)
+    std::map<art::Ptr<SimParticle>, std::pair<int, std::pair<float, float>>> info; // N(hits), (earliest time, p)
+    for (const auto& mc : sdmcs) {
+      const auto& sgs = mc.earlyStrawGasStep();
+      if (sgs.isNull()) continue;
+      const auto& simp = sgs->simParticle();
+      if (std::abs(simp->pdgId()) != 11 && std::abs(simp->pdgId()) != 13) continue; // e+- and mu+-
+      auto& in = info[simp];
+      if (in.first == 0 || sgs->time() < in.second.first) in.second = {float(sgs->time()), float(sgs->momvec().mag())};
+      ++in.first;
+    }
+    art::Ptr<SimParticle> sig;
+    _nsh = 0;
+    float pbest = 0.;
+    for (const auto& [simp, in] : info) {
+      if (in.first < _minHits || in.second.second <= pbest) continue;
+      pbest = in.second.second;
+      sig = simp;
+      _nsh = in.first;
+    }
+
+    // truth at the earliest step of the signal particle
+    _pdg = 0; _p = _pt = _pz = _rtrue = _lambdatrue = _xctrue = _yctrue = _d0true = 0.;
+    float tmin = 1.e12;
+    for (const auto& mc : sdmcs) {
+      const auto& sgs = mc.earlyStrawGasStep();
+      if (sgs.isNull() || sgs->simParticle() != sig || sgs->time() >= tmin) continue;
+      tmin = sgs->time();
+      const auto mom = sgs->momvec();
+      const auto pos = sgs->position();
+      _pdg = sig->pdgId();
+      _p = mom.mag(); _pt = mom.perp(); _pz = mom.z();
+      const float q = (sig->pdgId() == 11 || sig->pdgId() == 13 || sig->pdgId() == -211 || sig->pdgId() == -2212) ? -1. : 1.;
+      _rtrue = _pt / std::fabs(_bzConv);
+      _lambdatrue = _pz / std::fabs(_bzConv);
+      // circle center: pos + R * q * sign(B) * (py, -px)/pT
+      const float s = q * (_bzConv > 0 ? 1. : -1.);
+      _xctrue = pos.x() + _rtrue * s * mom.y() / _pt;
+      _yctrue = pos.y() - _rtrue * s * mom.x() / _pt;
+      _d0true = std::sqrt(_xctrue*_xctrue + _yctrue*_yctrue) - _rtrue;
+    }
+
+    // signal combo hits in the helix finder input collection: hits with a majority of signal straws
+    std::vector<float> zsig;
+    std::vector<size_t> sigHits; // signal combo hits in the input collection
+    std::vector<StrawHitIndex> shids;
+    for (size_t ich = 0; ich < chcol.size(); ++ich) {
+      shids.clear();
+      chcol.fillStrawHitIndices(ich, shids);
+      int nsig = 0;
+      for (auto sh : shids) {
+        const auto& sgs = sdmcs.at(sh).earlyStrawGasStep();
+        if (sgs.isNonnull() && sgs->simParticle() == sig) ++nsig;
+      }
+      if (2 * nsig > (int)shids.size()) { zsig.push_back(chcol[ich].pos().z()); sigHits.push_back(ich); }
+    }
+    _nchsig = zsig.size();
+    _sx.clear(); _sy.clear(); _sz.clear(); _snsh.clear(); _shel.clear();
+    for (auto ich : sigHits) {
+      const auto& hit = chcol[ich];
+      _sx.push_back(hit.pos().x()); _sy.push_back(hit.pos().y()); _sz.push_back(hit.pos().z());
+      _snsh.push_back(hit.nStrawHits());
+      int mask = 0;
+      for (size_t ih = 0; ih < hscol.size() && ih < 30; ++ih) {
+        for (const auto& hh : hscol[ih].hits()) {
+          if ((hh.pos() - hit.pos()).Mag2() < 1.e-6) { mask |= (1 << ih); break; }
+        }
+      }
+      _shel.push_back(mask);
+    }
+    _nsig3true = 0;
+    for (auto ich : sigHits) if (circleChi2(chcol[ich], _xctrue, _yctrue, _rtrue) < 9.f) ++_nsig3true;
+    zInfo(zsig, _zminsig, _zmaxsig, _gapsig, _nzsig);
+
+    // helices
+    _hlambda.clear(); _hradius.clear(); _hxc.clear(); _hyc.clear(); _hd0.clear(); _hzmin.clear(); _hzmax.clear();
+    _hgap.clear(); _hchi2xy.clear(); _hchi2zphi.clear(); _hnsh.clear(); _hnshsig.clear(); _hnz.clear();
+    _hchi2creco.clear(); _hchi2ctrue.clear(); _hnsig3reco.clear(); _hnsig3true.clear();
+    for (const auto& hs : hscol) {
+      const auto& h = hs.helix();
+      _hlambda.push_back(h.lambda());
+      _hradius.push_back(h.radius());
+      _hxc.push_back(h.centerx());
+      _hyc.push_back(h.centery());
+      _hd0.push_back(h.rcent() - h.radius());
+      _hchi2xy.push_back(h.chi2dXY());
+      _hchi2zphi.push_back(h.chi2dZPhi());
+      const auto& hits = hs.hits();
+      std::vector<float> z;
+      int nsh = 0, nshsig = 0;
+      for (size_t i = 0; i < hits.size(); ++i) {
+        z.push_back(hits[i].pos().z());
+        shids.clear();
+        hits.fillStrawHitIndices(i, shids);
+        for (auto sh : shids) {
+          ++nsh;
+          const auto& sgs = sdmcs.at(sh).earlyStrawGasStep();
+          if (sgs.isNonnull() && sgs->simParticle() == sig) ++nshsig;
+        }
+      }
+      float zmin, zmax, gap;
+      int nz;
+      zInfo(z, zmin, zmax, gap, nz);
+      _hzmin.push_back(zmin); _hzmax.push_back(zmax); _hgap.push_back(gap); _hnz.push_back(nz);
+      _hnsh.push_back(nsh); _hnshsig.push_back(nshsig);
+      float c2r = 0., c2t = 0.;
+      for (size_t i = 0; i < hits.size(); ++i) {
+        c2r += circleChi2(hits[i], h.centerx(), h.centery(), h.radius());
+        c2t += circleChi2(hits[i], _xctrue, _yctrue, _rtrue);
+      }
+      const float ndof = std::max<float>(1.f, hits.size() - 3.f);
+      _hchi2creco.push_back(c2r / ndof);
+      _hchi2ctrue.push_back(c2t / ndof);
+      int n3r = 0, n3t = 0;
+      for (auto ich : sigHits) {
+        if (circleChi2(chcol[ich], h.centerx(), h.centery(), h.radius()) < 9.f) ++n3r;
+        if (circleChi2(chcol[ich], _xctrue, _yctrue, _rtrue) < 9.f) ++n3t;
+      }
+      _hnsig3reco.push_back(n3r);
+      _hnsig3true.push_back(n3t);
+    }
+    _tree->Fill();
+  }
+
+} // namespace mu2e
+
+DEFINE_ART_MODULE(mu2e::AprHelixTruthAna)

From 8e8b3d1af957247013ed18a0d712ec009d99c094 Mon Sep 17 00:00:00 2001
From: michaelmackenzie <michaelmackenzie@u.northwestern.edu>
Date: Mon, 19 Oct 2026 10:01:15 -0500
Subject: [PATCH 4/5] HelixTool: out-of-panel crossings keep weight 1;
 ReadTriggerInfo useTrackerStatus

- liveFraction: a face crossing with no straw of that face at the point (radially
  outside the panel's straws, or beyond the straw ends, e.g. between the panels of
  the face) keeps the old weight 1 instead of taking the end straws' status.
- ReadTriggerInfo: add useTrackerStatus (default true, as in HelixFilter) so the
  hit-ratio histograms can match a filter run without the tracker status.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 Mu2eUtilities/inc/HelixTool.hh        |  3 ++-
 Mu2eUtilities/src/HelixTool.cc        |  9 +++++++++
 Trigger/src/ReadTriggerInfo_module.cc | 11 ++++++++---
 3 files changed, 19 insertions(+), 4 deletions(-)

diff --git a/Mu2eUtilities/inc/HelixTool.hh b/Mu2eUtilities/inc/HelixTool.hh
index 69582ab65d..370c5e78ca 100644
--- a/Mu2eUtilities/inc/HelixTool.hh
+++ b/Mu2eUtilities/inc/HelixTool.hh
@@ -44,7 +44,8 @@ namespace mu2e {
     //expected intersections of the helix with the tracker plane faces inside the tracker
     //active radius. If a StrawNoSignal function is given, each face crossing is weighted by
     //the fraction of live straws (nearest straw in each layer) at the crossing point, so that
-    //dead panels and straws don't lower the ratio. Evaluated on the first call.
+    //dead panels and straws don't lower the ratio. A crossing with no straw of the face at that
+    //point (outside the panel's straws) keeps the weight 1. Evaluated on the first call.
     float  hitRatio         () const;
 
     void   dirOfProp(float& slope, float& slopeErr, float& chi2ndof);
diff --git a/Mu2eUtilities/src/HelixTool.cc b/Mu2eUtilities/src/HelixTool.cc
index e33982dd25..afb9eb712b 100644
--- a/Mu2eUtilities/src/HelixTool.cc
+++ b/Mu2eUtilities/src/HelixTool.cc
@@ -142,8 +142,17 @@ namespace mu2e {
     const int nStraws = best->nStraws();
     float rho0   = strawRho(0);
     float rhoMax = strawRho(nStraws-1);
+    // no straw of this face at the crossing (radially outside the panel, or beyond the straw
+    // ends, e.g. between the panels of the face): keep the weight without the status
+    if (bestRho < rho0 || bestRho > rhoMax) return 1.;
     int   sEst   = (rhoMax > rho0) ? std::lround((bestRho-rho0)/(rhoMax-rho0)*(nStraws-1)) : 0;
     sEst = std::clamp(sEst, 0, nStraws-1);
+    {
+      const Straw& straw = best->getStraw(sEst);
+      const auto&  mid   = straw.getMidPoint();
+      const auto&  dir   = straw.getDirection();
+      if (std::fabs((X-mid.x())*dir.x() + (Y-mid.y())*dir.y()) > straw.halfLength()) return 1.;
+    }
 
     float nLive(0);
     for (int layer=0; layer<StrawId::_nlayers; ++layer) {
diff --git a/Trigger/src/ReadTriggerInfo_module.cc b/Trigger/src/ReadTriggerInfo_module.cc
index d69525bcec..9b673a47d3 100644
--- a/Trigger/src/ReadTriggerInfo_module.cc
+++ b/Trigger/src/ReadTriggerInfo_module.cc
@@ -194,6 +194,7 @@ namespace mu2e {
       fhicl::Atom<art::InputTag>     genCountTag        {Name("genCount"           ), Comment("GenEventCount label"                       ), "genCounter" };
       fhicl::Atom<art::InputTag>     PBITag             {Name("PBITag"             ), Comment("ProtonBunchIntensity label"                ), "PBISim" };
       fhicl::Atom<string>            processName        {Name("processName"        ), Comment("globalTrigger"                             ), ""  };
+      fhicl::Atom<bool>              useTrackerStatus   {Name("useTrackerStatus"   ), Comment("don't count straws without signal (TrackerStatus) as expected hits in the helix hit ratio; set as in the HelixFilter"), true };
     };
 
     explicit ReadTriggerInfo(const art::EDAnalyzer::Table<Config>& config);
@@ -263,7 +264,8 @@ namespace mu2e {
 
     const mu2e::Tracker*      _tracker = nullptr;
     ProditionsHandle<TrackerStatus> _trackerStatus_h;
-    HelixTool::StrawNoSignal  _noSignal;  // straws not counted as expected hits in the helix hit ratio
+    bool                      _useTrackerStatus = true;
+    HelixTool::StrawNoSignal  _noSignal;  // straws not counted as expected hits in the helix hit ratio; empty: all straws live
 
     float  _minPOT, _maxPOT;
     const int _genOccIndex;
@@ -288,6 +290,7 @@ namespace mu2e {
     _maxPOT              (3.5e8                      ),
     _genOccIndex         (2*_nTrackTrig + _nCaloTrig )
   {
+    _useTrackerStatus = config().useTrackerStatus();
 
     _trigAll.        resize(_nMaxTrig);
     _trigFinal.      resize(_nMaxTrig);
@@ -770,8 +773,10 @@ namespace mu2e {
   //--------------------------------------------------------------------------------//
   void ReadTriggerInfo::analyze(const art::Event& event) {
     // same definition of the helix hit ratio as in HelixFilter
-    const TrackerStatus* status = &_trackerStatus_h.get(event.id());
-    _noSignal = [status](StrawId const& sid) { return status->noSignal(sid) || status->suppress(sid); };
+    if (_useTrackerStatus) {
+      const TrackerStatus* status = &_trackerStatus_h.get(event.id());
+      _noSignal = [status](StrawId const& sid) { return status->noSignal(sid) || status->suppress(sid); };
+    }
 
     if(!_useNGen) _nProcess += 1;
 

From a336faff6d6e1571cf1b1a1f14415ac1def0e11a Mon Sep 17 00:00:00 2001
From: michaelmackenzie <michaelmackenzie@u.northwestern.edu>
Date: Mon, 19 Oct 2026 10:08:11 -0500
Subject: [PATCH 5/5] AgnosticHelixFinder: address review of the
 best-of-K/alias/target-refit options

- resolveGapAlias: refit the current line over the same points as the +-2pi
  alternatives, with a used calo cluster in all of them (as in recoverPoints).
- targetRefit/scanDphiDz: keep a used calo cluster in both the circle and the
  line fits, as in the other fits.
- A timeout during the triplet search saves the best stored candidate instead
  of dropping it.
- The K = 1 save path calls saveAndContinue, so the "search for another helix"
  rule has one home.
- The six new parameters are mandatory; their values (all options off) are set
  in the AgnosticHelixFinder table of CalPatRec/fcl/prolog.fcl.
- AprHelixTruthAna: PDGCode names, the _bzConv unit comment, and a note that the
  straw hit index is the StrawDigiMC index while FilterHits is false.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 Analyses/src/AprHelixTruthAna_module.cc     |  8 ++-
 CalPatRec/fcl/prolog.fcl                    |  6 ++
 CalPatRec/src/AgnosticHelixFinder_module.cc | 61 ++++++++++++---------
 3 files changed, 46 insertions(+), 29 deletions(-)

diff --git a/Analyses/src/AprHelixTruthAna_module.cc b/Analyses/src/AprHelixTruthAna_module.cc
index b07b36af77..a713f860aa 100644
--- a/Analyses/src/AprHelixTruthAna_module.cc
+++ b/Analyses/src/AprHelixTruthAna_module.cc
@@ -22,6 +22,7 @@
 #include "fhiclcpp/types/Atom.h"
 
 #include "Offline/BFieldGeom/inc/BFieldManager.hh"
+#include "Offline/DataProducts/inc/PDGCode.hh"
 #include "Offline/GeometryService/inc/DetectorSystem.hh"
 #include "Offline/GeometryService/inc/GeomHandle.hh"
 #include "Offline/MCDataProducts/inc/StrawDigiMC.hh"
@@ -66,7 +67,7 @@ namespace mu2e {
     int _minHits;
     std::vector<std::string> _paths;
     std::vector<int> _pass; // per entry of _paths
-    float _bzConv = 0.; // mm per MeV/c
+    float _bzConv = 0.; // MeV/c per mm
     TTree* _tree = nullptr;
 
     // event
@@ -207,7 +208,7 @@ namespace mu2e {
       const auto& sgs = mc.earlyStrawGasStep();
       if (sgs.isNull()) continue;
       const auto& simp = sgs->simParticle();
-      if (std::abs(simp->pdgId()) != 11 && std::abs(simp->pdgId()) != 13) continue; // e+- and mu+-
+      if (std::abs(simp->pdgId()) != PDGCode::e_minus && std::abs(simp->pdgId()) != PDGCode::mu_minus) continue; // e+- and mu+-
       auto& in = info[simp];
       if (in.first == 0 || sgs->time() < in.second.first) in.second = {float(sgs->time()), float(sgs->momvec().mag())};
       ++in.first;
@@ -233,7 +234,7 @@ namespace mu2e {
       const auto pos = sgs->position();
       _pdg = sig->pdgId();
       _p = mom.mag(); _pt = mom.perp(); _pz = mom.z();
-      const float q = (sig->pdgId() == 11 || sig->pdgId() == 13 || sig->pdgId() == -211 || sig->pdgId() == -2212) ? -1. : 1.;
+      const float q = (sig->pdgId() == PDGCode::e_minus || sig->pdgId() == PDGCode::mu_minus) ? -1. : 1.; // only e+- and mu+- reach here
       _rtrue = _pt / std::fabs(_bzConv);
       _lambdatrue = _pz / std::fabs(_bzConv);
       // circle center: pos + R * q * sign(B) * (py, -px)/pT
@@ -251,6 +252,7 @@ namespace mu2e {
       shids.clear();
       chcol.fillStrawHitIndices(ich, shids);
       int nsig = 0;
+      // the straw hit index is the StrawDigiMC index: true while the hit makers don't filter (FilterHits false)
       for (auto sh : shids) {
         const auto& sgs = sdmcs.at(sh).earlyStrawGasStep();
         if (sgs.isNonnull() && sgs->simParticle() == sig) ++nsig;
diff --git a/CalPatRec/fcl/prolog.fcl b/CalPatRec/fcl/prolog.fcl
index 4c21c55218..851b984ce5 100644
--- a/CalPatRec/fcl/prolog.fcl
+++ b/CalPatRec/fcl/prolog.fcl
@@ -546,6 +546,12 @@ CalPatRec : { @table::CalPatRec
             caloClusterSigma        : 30.0
             minNHelixComboHits      : 10
             minNHelixStrawHits      : 12
+            maxHelixCandidates      : 1    # save the first viable candidate
+            acceptHelixStrawHits    : 0    # (maxHelixCandidates > 1) save a candidate with this many straw hits at once; 0: never
+            minAliasGap             : 0.   # (mm) 2pi phi-z alias test across the largest z gap; 0: off
+            targetRefitSigma        : 0.   # (mm) target-constrained refit of the best candidate; 0: off
+            targetRefitPhiWindow    : 0.5  # (rad) target refit dphi/dz scan window
+            targetRefitDphiDzBins   : 100  # target refit dphi/dz scan bins per helicity
             minHelixPerpMomentum    : 20.0
             maxHelixPerpMomentum    : 140.0
             minHelixMomentum        : 30.0
diff --git a/CalPatRec/src/AgnosticHelixFinder_module.cc b/CalPatRec/src/AgnosticHelixFinder_module.cc
index f51eaaa653..6b82bf3e91 100644
--- a/CalPatRec/src/AgnosticHelixFinder_module.cc
+++ b/CalPatRec/src/AgnosticHelixFinder_module.cc
@@ -94,12 +94,12 @@ namespace mu2e {
       fhicl::Atom<float>           caloClusterSigma       {Name("caloClusterSigma"     ), Comment("error assigned to calo clust")  };
       fhicl::Atom<int>             minNHelixStrawHits     {Name("minNHelixStrawHits"   ), Comment("straw hit save threshold"    )  };
       fhicl::Atom<int>             minNHelixComboHits     {Name("minNHelixComboHits"   ), Comment("combo hit save threshold"    )  };
-      fhicl::Atom<int>             maxHelixCandidates     {Name("maxHelixCandidates"   ), Comment("save the best (most straw hits) of up to this many candidates per search; 1: save the first"), 1};
-      fhicl::Atom<int>             acceptHelixStrawHits   {Name("acceptHelixStrawHits" ), Comment("save a candidate with at least this many straw hits without searching further; <= 0: never"), 0};
-      fhicl::Atom<float>           minAliasGap            {Name("minAliasGap"          ), Comment("test the 2pi phi-z alias across the largest z gap of a helix's hits if it is at least this (mm); <= 0: off"), 0.f};
-      fhicl::Atom<float>           targetRefitSigma       {Name("targetRefitSigma"     ), Comment("if the best of the maxHelixCandidates (> 1) candidates has |d0| > this (mm), refit its circle through the stopping target (error = this) with a dphi/dz scan; the refit replaces it if it has at least as many straw hits; <= 0: off"), 0.f};
-      fhicl::Atom<float>           targetRefitPhiWindow   {Name("targetRefitPhiWindow" ), Comment("phi window (rad) for hits to count on a line in the target refit's dphi/dz scan"), 0.5f};
-      fhicl::Atom<int>             targetRefitDphiDzBins  {Name("targetRefitDphiDzBins"), Comment("dphi/dz bins per helicity in the target refit's dphi/dz scan"), 100};
+      fhicl::Atom<int>             maxHelixCandidates     {Name("maxHelixCandidates"   ), Comment("save the best (most straw hits) of up to this many candidates per search; 1: save the first")};
+      fhicl::Atom<int>             acceptHelixStrawHits   {Name("acceptHelixStrawHits" ), Comment("save a candidate with at least this many straw hits without searching further; <= 0: never")};
+      fhicl::Atom<float>           minAliasGap            {Name("minAliasGap"          ), Comment("test the 2pi phi-z alias across the largest z gap of a helix's hits if it is at least this (mm); <= 0: off")};
+      fhicl::Atom<float>           targetRefitSigma       {Name("targetRefitSigma"     ), Comment("if the best of the maxHelixCandidates (> 1) candidates has |d0| > this (mm), refit its circle through the stopping target (error = this) with a dphi/dz scan; the refit replaces it if it has at least as many straw hits; <= 0: off")};
+      fhicl::Atom<float>           targetRefitPhiWindow   {Name("targetRefitPhiWindow" ), Comment("phi window (rad) for hits to count on a line in the target refit's dphi/dz scan")};
+      fhicl::Atom<int>             targetRefitDphiDzBins  {Name("targetRefitDphiDzBins"), Comment("dphi/dz bins per helicity in the target refit's dphi/dz scan")};
       fhicl::Atom<float>           minHelixPerpMomentum   {Name("minHelixPerpMomentum" ), Comment("min pt of helix"             )  };
       fhicl::Atom<float>           maxHelixPerpMomentum   {Name("maxHelixPerpMomentum" ), Comment("max pt of helix"             )  };
       fhicl::Atom<float>           minHelixMomentum       {Name("minHelixMomentum"     ), Comment("min momentum of helix"       )  };
@@ -866,7 +866,8 @@ namespace mu2e {
       if(shouldExitForTimeout(__func__)) {
         if(_doTiming > 0) _watch->StopTime(__func__);
         if(_doTiming > 1) _watch->StopTime("triplet-i");
-        return false;
+        // keep the best candidate found before the timeout (maxHelixCandidates > 1)
+        return _candidates.empty() ? false : saveCandidates(tc, HSColl);
       }
       bool uselessSeed = true;
       setTripletI(i, tripletInfo, loopCondition);
@@ -946,14 +947,8 @@ namespace mu2e {
             if(_doTiming > 3) _watch->StopTime("triplet-k");
             return saveCandidates(tc, HSColl);
           } else if (nStrawHitsInHelix >= _minNHelixStrawHits && nComboHitsInHelix >= _minNHelixComboHits) {
-            saveHelix(tc, HSColl);
-            if (_diagLevel > 0) { _diagInfo.nHelices++; }
-            // we only want to search for another helix if we have enough remaining hits after saving helix
-            const int remainingStrawHits = nStrawHitsInTimeCluster - nStrawHitsInHelix;
-            const int remainingComboHits = nComboHitsInTimeCluster - nComboHitsInHelix;
-            const bool findAnotherHelix = remainingStrawHits >= _minNHelixStrawHits && remainingComboHits >= _minNHelixComboHits;
-            if(findAnotherHelix && _doAverageFlag) resetFlags(); // reset flags if needed
-
+            const bool findAnotherHelix = saveAndContinue(tc, HSColl, nStrawHitsInHelix, nComboHitsInHelix,
+                                                          nStrawHitsInTimeCluster, nComboHitsInTimeCluster);
             // Helix found, so we can exit this search
             if(_doTiming > 0) _watch->StopTime(__func__);
             if(_doTiming > 1) _watch->StopTime("triplet-i");
@@ -1000,7 +995,7 @@ namespace mu2e {
       _circleFitter.addPoint(0.f, 0.f, wT);
       for (size_t i = 0; i < _tcHits.size(); i++) {
         const auto& h = _tcHits[i];
-        if (!h.used || h.inHelix || h.hitIndice == HitType::STOPPINGTARGET || h.hitIndice == HitType::CALOCLUSTER) continue;
+        if (!h.used || h.inHelix || h.hitIndice == HitType::STOPPINGTARGET) continue; // a used calo cluster is kept, as in the other fits
         computeCircleError2(i, xC, yC);
         _circleFitter.addPoint(getPos(i).x(), getPos(i).y(), 1.f / _tcHits[i].circleError2);
       }
@@ -1062,6 +1057,7 @@ namespace mu2e {
                                             int nStrawHitsInTimeCluster, int nComboHitsInTimeCluster) {
     saveHelix(tc, HSColl);
     if (_diagLevel > 0) { _diagInfo.nHelices++; }
+    // we only want to search for another helix if we have enough remaining hits after saving helix
     const int remainingStrawHits = nStrawHitsInTimeCluster - nStrawHitsInHelix;
     const int remainingComboHits = nComboHitsInTimeCluster - nComboHitsInHelix;
     const bool findAnotherHelix = remainingStrawHits >= _minNHelixStrawHits && remainingComboHits >= _minNHelixComboHits;
@@ -1725,12 +1721,14 @@ namespace mu2e {
   //-----------------------------------------------------------------------------
   bool AgnosticHelixFinder::resolveGapAlias() {
     if(_doTiming > 4) _watch->SetTime(__func__);
-    // used tracker hits, in decreasing z (the _tcHits order)
-    std::vector<size_t> idx;
+    // used tracker hits, in decreasing z (the _tcHits order); a used calo cluster is in the line
+    // fits, as in recoverPoints, on the high-z side of any tracker gap, but not in the gap search
+    std::vector<size_t> idx, calo;
     for (size_t i = 0; i < _tcHits.size(); i++) {
       const auto& h = _tcHits[i];
-      if (!h.used || h.inHelix || h.hitIndice == HitType::STOPPINGTARGET || h.hitIndice == HitType::CALOCLUSTER) continue;
-      idx.push_back(i);
+      if (!h.used || h.inHelix || h.hitIndice == HitType::STOPPINGTARGET) continue;
+      if (h.hitIndice == HitType::CALOCLUSTER) calo.push_back(i);
+      else                                     idx.push_back(i);
     }
     if (idx.size() < 4) { if(_doTiming > 4) _watch->StopTime(__func__); return false; }
     size_t igap = 0; // hits idx[0..igap] are on the high-z side
@@ -1741,17 +1739,23 @@ namespace mu2e {
     }
     if (gap < _minAliasGap) { if(_doTiming > 4) _watch->StopTime(__func__); return false; }
 
+    // the current line (shift 0) is refit over the same points as the alternatives
     const float radius = _circleFitter.radius();
-    float bestChi2 = _lineFitter.chi2Dof();
+    float bestChi2 = 0.f;
     int   bestShift = 0;
     ::LsqSums2 bestFitter;
-    for (int shift : {-1, 1}) {
+    for (int shift : {0, -1, 1}) {
       ::LsqSums2 fitter;
       for (size_t k = 0; k < idx.size(); ++k) {
         const auto& h = _tcHits[idx[k]];
         const int n = h.helixPhiCorrection + (k > igap ? shift : 0);
         fitter.addPoint(getPos(idx[k]).z(), h.helixPhi + _twopi * n, 1.f / h.helixPhiError2);
       }
+      for (auto i : calo) {
+        const auto& h = _tcHits[i];
+        fitter.addPoint(getPos(i).z(), h.helixPhi + _twopi * h.helixPhiCorrection, 1.f / h.helixPhiError2);
+      }
+      if (shift == 0) { bestChi2 = fitter.chi2Dof(); continue; }
       const float mom2 = computeHelixMomentum2(radius, fitter.dydx());
       if (mom2 < _minHelixMomentumSq || mom2 > _maxHelixMomentumSq) continue;
       if (fitter.chi2Dof() < bestChi2) {
@@ -1847,11 +1851,16 @@ namespace mu2e {
       _lineFitter.addPoint(z[k], phi[k] + _twopi * n, 1.f / h.helixPhiError2);
       _circleFitter.addPoint(getPos(cand[k]).x(), getPos(cand[k]).y(), 1.f / h.circleError2);
     }
-    // keep the calo cluster in the circle fit, as in the other fits
+    // keep a used calo cluster in the circle and line fits, as in the other fits
     for (size_t i = 0; i < _tcHits.size(); i++) {
-      const auto& h = _tcHits[i];
-      if (h.hitIndice == HitType::CALOCLUSTER && h.used && !h.inHelix)
-        _circleFitter.addPoint(getPos(i).x(), getPos(i).y(), 1.f / h.circleError2);
+      auto& h = _tcHits[i];
+      if (h.hitIndice != HitType::CALOCLUSTER || !h.used || h.inHelix) continue;
+      computeCircleError2(i, xC, yC);
+      computeHelixPhi(i, xC, yC);
+      const float zC = getPos(i).z();
+      h.helixPhiCorrection = std::round((bestS * zC + bestPhi0 - h.helixPhi) / _twopi);
+      _circleFitter.addPoint(getPos(i).x(), getPos(i).y(), 1.f / h.circleError2);
+      _lineFitter.addPoint(zC, h.helixPhi + _twopi * h.helixPhiCorrection, 1.f / h.helixPhiError2);
     }
     return done(true);
   }