diff --git a/Tr/TrackMonitors/CMakeLists.txt b/Tr/TrackMonitors/CMakeLists.txt index fffb741970034c08d70216e182d986be827493e4..609ad38ea7c0b99004fc59d6bb16f6b639785109 100644 --- a/Tr/TrackMonitors/CMakeLists.txt +++ b/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 diff --git a/Tr/TrackMonitors/src/UTHitEfficiencyMonitor.cpp b/Tr/TrackMonitors/src/UTHitEfficiencyMonitor.cpp new file mode 100644 index 0000000000000000000000000000000000000000..9a86f74a8de9429f635f6d80faf243079894ad2b --- /dev/null +++ b/Tr/TrackMonitors/src/UTHitEfficiencyMonitor.cpp @@ -0,0 +1,234 @@ +/*****************************************************************************\ +* (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 + +namespace LHCb { + + class UTHitEfficiencyMonitor + : public LHCb::Algorithm::Consumer> { + + private: + Gaudi::Property m_layerUnderStudy{ this, "LayerUnderStudy" }; + Gaudi::Property m_maxDoca{ this, "MaxDoca", 1 * Gaudi::Units::mm }; + Gaudi::Property m_maxTrackCov{ this, "MaxTrackCov", 0.5 * Gaudi::Units::mm }; + Gaudi::Property m_trackMaxChi2PerDoF{ this, "MaxTrackChi2PerDof", 3. }; + Gaudi::Property m_trackMinP{ this, "MinTrackP", 2000 * Gaudi::Units::MeV }; + Gaudi::Property m_trackMinPt{ this, "MinTrackPT", 400 * Gaudi::Units::MeV }; + + ToolHandle 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 m_layerUnderStudy_error{ this, "Failed to find UT layer " }; + mutable Gaudi::Accumulators::MsgCounter m_interpolator_warning{ this, "Failed to interpolate track" }; + mutable Gaudi::Accumulators::MsgCounter m_covariance_debug{ this, "Covariance too large" }; + mutable Gaudi::Accumulators::MsgCounter 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 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().cast(); + const float hit_yb = hit.get().cast(); + const float hit_ye = hit.get().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( hitpos.beginPoint(), hitpos.endPoint() ); + + auto trackLine = + Gaudi::Math::Line( 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