Merge branch 'test_UT_eff_FW' into 'master'

UT hit efficiency monitor

See merge request !3925

Tr/TrackMonitors/CMakeLists.txt

@@ -34,7 +34,8 @@ gaudi_add_module(TrackMonitors
         src/TrackTune.cpp
         src/TrackVertexMonitor.cpp
         src/TrackVPOverlapMonitor.cpp
-        src/UTTrackMonitor.cpp
+        src/UTHitEfficiencyMonitor.cpp
+	src/UTTrackMonitor.cpp
         src/VPTrackMonitor.cpp
         src/VPHitEfficiencyMonitor.cpp
         src/VertexCompare.cpp

Tr/TrackMonitors/src/UTHitEfficiencyMonitor.cpp

+/*****************************************************************************\
+* (c) Copyright 2024 CERN for the benefit of the LHCb Collaboration      *
+*                                                                             *
+* This software is distributed under the terms of the GNU General Public      *
+* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   *
+*                                                                             *
+* In applying this licence, CERN does not waive the privileges and immunities *
+* granted to it by virtue of its status as an Intergovernmental Organization  *
+* or submit itself to any jurisdiction.                                       *
+\*****************************************************************************/
+#include "LHCbAlgs/Consumer.h"
+#include "LHCbAlgs/Traits.h"
+
+#include "DetDesc/DetectorElement.h"
+#include "Detector/UT/ChannelID.h"
+#include "Detector/UT/UTConstants.h"
+#include "Event/PrHits.h"
+
+#include "UTDet/DeUTDetector.h"
+#include "UTDet/DeUTFace.h"
+#include "UTDet/DeUTLayer.h"
+#include "UTDet/DeUTModule.h"
+#include "UTDet/DeUTSector.h"
+#include "UTDet/DeUTSensor.h"
+#include "UTDet/DeUTSide.h"
+#include "UTDet/DeUTStation.h"
+#include "UTDet/DeUTStave.h"
+
+#include "GaudiKernel/SystemOfUnits.h"
+
+#include "Event/PrHits.h"
+#include "Event/Track.h"
+
+#include "Gaudi/Accumulators/StaticHistogram.h"
+
+#include "TrackInterfaces/ITrackInterpolator.h"
+
+#include "GaudiKernel/Vector3DTypes.h"
+#include "LHCbMath/GeomFun.h"
+#include "LHCbMath/Line.h"
+#include <cmath>
+
+namespace LHCb {
+
+  class UTHitEfficiencyMonitor
+      : public LHCb::Algorithm::Consumer<void( const LHCb::Track::Range&, const LHCb::Pr::UT::Hits&,
+                                               const DetectorElement&, const DeUTDetector& ),
+                                         LHCb::Algorithm::Traits::usesConditions<DetectorElement, DeUTDetector>> {
+
+  private:
+    Gaudi::Property<unsigned int> m_layerUnderStudy{ this, "LayerUnderStudy" };
+    Gaudi::Property<float>        m_maxDoca{ this, "MaxDoca", 1 * Gaudi::Units::mm };
+    Gaudi::Property<float>        m_maxTrackCov{ this, "MaxTrackCov", 0.5 * Gaudi::Units::mm };
+    Gaudi::Property<float>        m_trackMaxChi2PerDoF{ this, "MaxTrackChi2PerDof", 3. };
+    Gaudi::Property<float>        m_trackMinP{ this, "MinTrackP", 2000 * Gaudi::Units::MeV };
+    Gaudi::Property<float>        m_trackMinPt{ this, "MinTrackPT", 400 * Gaudi::Units::MeV };
+
+    ToolHandle<ITrackInterpolator> m_interpolator{ this, "Interpolator", "TrackInterpolator" };
+
+    mutable Gaudi::Accumulators::StaticHistogram<1> m_layer_Hit_recom{
+        this, "ReconstructionHit", "Reconstruction;HitPerLayer;Entries", { 10, -0.5, 9.5 } };
+
+    mutable Gaudi::Accumulators::StaticHistogram<1> m_matchedHitsPerTrack{
+        this, "MatchedHitsPerTrack", "Matched hits per Track;Hits;Entries", { 10, -0.5, 9.5 } };
+
+    mutable Gaudi::Accumulators::StaticProfileHistogram<2> m_efficiencyVsXY{
+        this,
+        "EfficiencyVsXY",
+        "Efficiency vs XY;X [mm];Y [mm]",
+        { 76, -908.2 * Gaudi::Units::mm, 908.2 * Gaudi::Units::mm },
+        { 32, -764.8 * Gaudi::Units::mm, 764.8 * Gaudi::Units::mm } };
+
+    mutable Gaudi::Accumulators::StaticHistogram<2> m_hitXY{
+        this,
+        "hitXY",
+        "XY;X [mm];Y [mm]",
+        { 76, -908.2 * Gaudi::Units::mm, 908.2 * Gaudi::Units::mm },
+        { 32, -764.8 * Gaudi::Units::mm, 764.8 * Gaudi::Units::mm } };
+    mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_efficiencyVsX{
+        this,
+        "EfficiencyVsX",
+        "Efficiency vs X;X [mm];Efficiency",
+        { 76, -908.2 * Gaudi::Units::mm, 908.2 * Gaudi::Units::mm } };
+
+    mutable Gaudi::Accumulators::StaticHistogram<2> m_All{
+        this,
+        "All Track",
+        "All Track;X [mm];Y [mm]",
+        { 76, -908.2 * Gaudi::Units::mm, 908.2 * Gaudi::Units::mm },
+        { 32, -764.8 * Gaudi::Units::mm, 764.8 * Gaudi::Units::mm } };
+
+    mutable Gaudi::Accumulators::StaticHistogram<2> m_covarianceVsXY{
+        this,
+        "CovarianceVsXY",
+        "Covariance vs XY;X [mm];Y [mm];Hits",
+        { 76, -908.2 * Gaudi::Units::mm, 908.2 * Gaudi::Units::mm },
+        { 32, -764.8 * Gaudi::Units::mm, 764.8 * Gaudi::Units::mm } };
+
+    mutable Gaudi::Accumulators::StaticHistogram<1> m_doca{
+        this, "Doca", "DOCA;DOCA [mm];Entries", { 100, -2.0 * Gaudi::Units::mm, 2.0 * Gaudi::Units::mm } };
+
+    mutable Gaudi::Accumulators::StaticHistogram<1> m_docaX{
+        this, "DocaX", "DOCA X;DOCA X [mm];Entries", { 100, -2.0 * Gaudi::Units::mm, 2.0 * Gaudi::Units::mm } };
+
+    mutable Gaudi::Accumulators::StaticHistogram<1> m_docaY{
+        this, "DocaY", "DOCA Y;DOCA Y [mm];Entries", { 100, -2.0 * Gaudi::Units::mm, 2.0 * Gaudi::Units::mm } };
+
+    mutable Gaudi::Accumulators::MsgCounter<MSG::ERROR>   m_layerUnderStudy_error{ this, "Failed to find UT layer " };
+    mutable Gaudi::Accumulators::MsgCounter<MSG::WARNING> m_interpolator_warning{ this, "Failed to interpolate track" };
+    mutable Gaudi::Accumulators::MsgCounter<MSG::DEBUG>   m_covariance_debug{ this, "Covariance too large" };
+    mutable Gaudi::Accumulators::MsgCounter<MSG::WARNING> track_hit_parallel{ this, "Track and hit are parallel" };
+
+  public:
+    UTHitEfficiencyMonitor( const std::string& name, ISvcLocator* pSvcLocator )
+        : Consumer( name, pSvcLocator,
+                    { KeyValue{ "TrackLocation", "" }, KeyValue{ "PrUTHitsLocation", "" },
+                      KeyValue{ "LHCbGeoLocation", LHCb::standard_geometry_top },
+                      KeyValue{ "UTDetectorLocation", DeUTDetLocation::location() } } ) {}
+
+    std::pair<double, double> correct_position_for_layer( int layer, double xt, double yt ) const {
+      if ( layer != 1 and layer != 2 ) return { xt, yt }; // not stereo
+
+      constexpr double sin_ster        = 0.087155743; // sin(5*TMath::Pi()/180.);
+      constexpr double cos_ster        = 0.99619470;  // cos(5*TMath::Pi()/180.);
+      constexpr double UT_sensor_width = 95.6 * Gaudi::Units::mm;
+
+      int sign = ( layer == 1 ? +1 : -1 ); // distinguish U and V
+
+      double xtc = xt - yt * ( sign * sin_ster ); // approximation
+      int    _si = ( xtc > 0.0 ) ? 1 : 0;
+      double axf = fabs( xtc ) / UT_sensor_width;
+      int    _st = int( axf );
+      double xst = ( ( _st + 0.5 ) * ( 2 * _si - 1 ) ) * UT_sensor_width;
+      double ytc = yt * cos_ster + ( xt - xst ) * ( sign * sin_ster );
+
+      return { xtc, ytc };
+    }
+
+    void operator()( const LHCb::Track::Range& tracks, const LHCb::Pr::UT::Hits& hits, const DetectorElement& lhcb,
+                     const DeUTDetector& utDet ) const override {
+
+      const auto& utlayergeom = utDet.getLayerGeom( m_layerUnderStudy );
+
+      if ( !utlayergeom.z ) {
+        ++m_layerUnderStudy_error;
+        return;
+      }
+
+      for ( const auto& track : tracks ) {
+
+        LHCb::State state;
+        if ( m_trackMaxChi2PerDoF.value() > 0 && track->chi2PerDoF() > m_trackMaxChi2PerDoF.value() ) continue;
+        if ( m_trackMinP.value() > 0 && track->p() < m_trackMinP.value() ) continue;
+        if ( m_trackMinPt.value() > 0 && track->pt() < m_trackMinPt.value() ) continue;
+
+        if ( !m_interpolator->interpolate( *track, utlayergeom.z, state, lhcb ).isSuccess() ) {
+          ++m_interpolator_warning;
+          continue;
+        }
+
+        m_covarianceVsXY[{ -state.x(), state.y() }] +=
+            ( std::sqrt( state.covariance()( 0, 0 ) ) > m_maxTrackCov.value() ) ? 1 : 0;
+
+        if ( std::sqrt( state.covariance()( 0, 0 ) ) > m_maxTrackCov.value() ) {
+          ++m_covariance_debug;
+          continue;
+        }
+
+        int   nHits    = 0;
+        float matchHit = 0;
+
+        for ( unsigned int index = 0; index < hits.size(); index++ ) {
+
+          if ( hits.id( index ).layer() != m_layerUnderStudy ) continue;
+
+          const auto hit = hits.scalar()[index];
+
+          const float hit_z  = hit.get<LHCb::Pr::UT::UTHitsTag::zAtYEq0>().cast();
+          const float hit_yb = hit.get<LHCb::Pr::UT::UTHitsTag::yBegin>().cast();
+          const float hit_ye = hit.get<LHCb::Pr::UT::UTHitsTag::yEnd>().cast();
+
+          const auto& channelid = hits.id( index );
+          const auto  hitpos    = utDet.trajectory( channelid, 0. );
+
+          LHCb::State temp_state;
+
+          if ( !m_interpolator->interpolate( *track, hit_z, temp_state, lhcb ).isSuccess() ) {
+            ++m_interpolator_warning;
+            continue;
+          }
+
+          auto hitLine = Gaudi::Math::Line<Gaudi::XYZPoint, Gaudi::XYZVector>( hitpos.beginPoint(), hitpos.endPoint() );
+
+          auto trackLine =
+              Gaudi::Math::Line<Gaudi::XYZPoint, Gaudi::XYZVector>( temp_state.position(), temp_state.slopes() );
+
+          Gaudi::XYZPoint trackPoint;
+          Gaudi::XYZPoint hitPoint;
+
+          if ( !Gaudi::Math::closestPoints( trackLine, hitLine, trackPoint, hitPoint ) ) {
+            ++track_hit_parallel;
+            continue;
+          }
+
+          if ( hitPoint.y() > std::max( hit_yb, hit_ye ) || hitPoint.y() < std::min( hit_yb, hit_ye ) ) continue;
+
+          auto docaVector = ( trackPoint - hitPoint );
+          auto signedDoca = std::copysign( docaVector.R(), docaVector.x() );
+
+          if ( std::abs( signedDoca ) <= m_maxDoca ) {
+            nHits++;
+            matchHit = 1;
+          }
+
+          ++m_doca[signedDoca];
+          ++m_docaX[docaVector.x()];
+          ++m_docaY[docaVector.y()];
+          ++m_hitXY[{ hitPoint.x(), hitPoint.y() }];
+        }
+
+        auto [xtc, ytc] = correct_position_for_layer( m_layerUnderStudy, state.x(), state.y() );
+
+        ++m_matchedHitsPerTrack[nHits];
+        ++m_layer_Hit_recom[track->nUTHits()];
+        m_efficiencyVsXY[{ xtc, ytc }] += matchHit;
+        m_efficiencyVsX[xtc] += matchHit;
+        m_All[{ xtc, ytc }] += 1;
+      }
+    }
+  };
+
+  DECLARE_COMPONENT_WITH_ID( UTHitEfficiencyMonitor, "UTHitEfficiencyMonitor" )
+
+} // namespace LHCb