HelixTool: don't count straws without signal in the helix hit ratio

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,44 @@ 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. 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);
 
   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

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,106 @@ 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);
+    // 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) {
+      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;
   }
 
 

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
 )

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"
@@ -192,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") };
     };
 
     explicit ReadTriggerInfo(const art::EDAnalyzer::Table<Config>& config);
@@ -260,6 +263,9 @@ namespace mu2e {
     occupancyHist_            _occupancyHist;
 
     const mu2e::Tracker*      _tracker = nullptr;
+    ProditionsHandle<TrackerStatus> _trackerStatus_h;
+    bool                      _useTrackerStatus;
+    HelixTool::StrawNoSignal  _noSignal;  // straws not counted as expected hits in the helix hit ratio; empty: all straws live
 
     float  _minPOT, _maxPOT;
     const int _genOccIndex;
@@ -284,6 +290,7 @@ namespace mu2e {
     _maxPOT              (3.5e8                      ),
     _genOccIndex         (2*_nTrackTrig + _nCaloTrig )
   {
+    _useTrackerStatus = config().useTrackerStatus();
 
     _trigAll.        resize(_nMaxTrig);
     _trigFinal.      resize(_nMaxTrig);
@@ -765,6 +772,11 @@ namespace mu2e {
 
   //--------------------------------------------------------------------------------//
   void ReadTriggerInfo::analyze(const art::Event& event) {
+    // same definition of the helix hit ratio as in HelixFilter
+    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;
 
@@ -1011,7 +1023,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);

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',

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
 )