diff --git a/Muon/MuonMonitors/src/MuonChamberMonitor.cpp b/Muon/MuonMonitors/src/MuonChamberMonitor.cpp index d7bf62fb4a203454d4f3cc2fcdaeb24ab193a4b5..eede85e6f8931ddca6d08445c95eec4481f62d1e 100644 --- a/Muon/MuonMonitors/src/MuonChamberMonitor.cpp +++ b/Muon/MuonMonitors/src/MuonChamberMonitor.cpp @@ -169,18 +169,18 @@ namespace LHCb { const MuonHitContainer& muonhits ) const override; private: - void fillHistograms( const DeMuonDetector& det, const DetectorElement& lhcb, const GeomCache& geometryinfo, - const LHCb::Particle::Range& probes, const LHCb::Particle::Range& tags, - const MuonHitContainer& muonhits ) const; - // The track extrapolator ToolHandle m_extrapolator = { this, "ReferenceExtrapolator", "TrackMasterExtrapolator" }; - std::tuple computeChi2( const GeomCache& geometryinfo, const xy_t& ext, - const xy_t& extMS, - const CommonMuonHitRange& hitContainer ) const; + struct chi2Result_t { + CommonMuonHit hit; + double chi2X; + double chi2Y; + }; + std::optional computeChi2( const GeomCache& geometryinfo, const xy_t& ext, const xy_t& extMS2, + const CommonMuonHitRange& hitContainer ) const; - int HitsInChi2( const GeomCache& geometryinfo, const xy_t& ext, const xy_t& extMS, + int HitsInChi2( const GeomCache& geometryinfo, const xy_t& ext, const xy_t& extMS2, const CommonMuonHitRange& hitContainer, int sign ) const; int xy2Chamber( const GeomCache& geometryinfo, const double x, const double y, unsigned int station, @@ -201,9 +201,12 @@ namespace LHCb { const LHCb::Track* m_track = nullptr; const LHCb::State* m_state = nullptr; - double m_probeP{ 0. }; - std::array m_probeChi2StX{ { 0., 0., 0., 0. } }; - std::array m_probeChi2StY{ { 0., 0., 0., 0. } }; + double m_probeP{ 0. }; + std::array m_probeChi2StX{ { 0., 0., 0., 0. } }; + std::array m_probeChi2StY{ { 0., 0., 0., 0. } }; + std::array m_constantsMS{ + { 220000., 377000., 565000., 785000. } }; // linear (X or Y) MS deviations times track momentum obtained from + // MC ArrayPairSt m_trackExt{ { { 0., 0. }, { 0., 0. }, { 0., 0. }, { 0., 0. } } }; ArrayPairSt m_trackMS{ { { 0., 0. }, { 0., 0. }, { 0., 0. }, { 0., 0. } } }; ArrayPairSt m_chit_XY{ { { 0., 0. }, { 0., 0. }, { 0., 0. }, { 0., 0. } } }; @@ -436,12 +439,6 @@ namespace LHCb { void MuonChamberMonitor::operator()( const DeMuonDetector& det, const DetectorElement& lhcb, const GeomCache& geometryinfo, const LHCb::Particle::Range& probes, const LHCb::Particle::Range& tags, const MuonHitContainer& muonhits ) const { - fillHistograms( det, lhcb, geometryinfo, probes, tags, muonhits ); - } - - void MuonChamberMonitor::fillHistograms( const DeMuonDetector& det, const DetectorElement& lhcb, - const GeomCache& geometryinfo, const LHCb::Particle::Range& probes, - const LHCb::Particle::Range& tags, const MuonHitContainer& muonhits ) const { auto daqHelper = det.getUpgradeDAQInfo(); auto& geometry = *lhcb.geometry(); unsigned int ntags = 0; @@ -464,8 +461,7 @@ namespace LHCb { candMuon.m_proto = probe->proto(); candMuon.m_track = probe->proto()->track(); candMuon.m_state = probe->proto()->track()->stateAt( LHCb::State::Location::LastMeasurement ); - // if (!candMuon.m_state) { continue; } - State trackState = closestState( *probe->proto()->track(), candMuon.m_state->z() ); + State trackState = closestState( *probe->proto()->track(), candMuon.m_state->z() ); // Match region of extrapolation for ( unsigned int s = 0; s != nMuonStations; ++s ) { StatusCode sc = m_extrapolator->propagate( trackState, geometryinfo.m_stationZ[s], geometry ); @@ -473,6 +469,8 @@ namespace LHCb { warning() << fmt::format( "Something went wrong in track extrapolation in station M{}", s + 2 ) << endmsg; candMuon.m_trackExt[s] = { trackState.x(), trackState.y() }; candMuon.m_trackMS[s] = { trackState.errX2(), trackState.errY2() }; + auto ms_error_sqared = pow( candMuon.m_constantsMS[s] / candMuon.m_probeP, 2 ); + candMuon.m_trackMS[s] = { ms_error_sqared, ms_error_sqared }; int region = 0; while ( region < 3 ) { if ( fabs( candMuon.m_trackExt[s].first ) < ( geometryinfo.m_regionOuterX[s][region] ) && @@ -483,23 +481,27 @@ namespace LHCb { candMuon.m_extrapReg[s] = region; candMuon.m_idReg[s] = s * nMuonStations + region; } + bool outOfAcceptance = false; for ( unsigned int s = 0; s != nMuonStations; ++s ) { - ++m_histos->m_probeExtrap[{ candMuon.m_trackExt[s].first, candMuon.m_trackExt[s].second }]; - auto [chit, chi2X, chi2Y] = + auto chi2_result = computeChi2( geometryinfo, candMuon.m_trackExt[s], candMuon.m_trackMS[s], muonhits.station( s ).hits() ); + if ( !chi2_result ) { + outOfAcceptance = true; + break; + } - candMuon.m_chitReg[s] = chit.region(); - candMuon.m_chit[s] = chit; - candMuon.m_chit_XY[s] = { chit.x(), chit.y() }; - candMuon.m_probeChi2StX[s] = chi2X; - candMuon.m_probeChi2StY[s] = chi2Y; - ++m_histos->m_Chi2StationX[s][candMuon.m_probeChi2StX[s]]; - ++m_histos->m_Chi2StationY[s][candMuon.m_probeChi2StY[s]]; + candMuon.m_chitReg[s] = chi2_result->hit.region(); + candMuon.m_chit[s] = chi2_result->hit; + candMuon.m_chit_XY[s] = { chi2_result->hit.x(), chi2_result->hit.y() }; + candMuon.m_probeChi2StX[s] = chi2_result->chi2X; + candMuon.m_probeChi2StY[s] = chi2_result->chi2Y; + ++m_histos->m_Chi2StationX[s][candMuon.m_probeChi2StX[s]]; // Chi2 histograms are filled for all the hits + ++m_histos->m_Chi2StationY[s][candMuon.m_probeChi2StY[s]]; // within acceptance, for debug purposes. }; + if ( outOfAcceptance ) continue; for ( unsigned int s = 0; s != nMuonStations; ++s ) { - bool otherStationsMatched = true; for ( unsigned ss = 0; ss != nMuonStations; ++ss ) { if ( ss == s ) continue; // Ignore the station we are looking at @@ -574,43 +576,51 @@ namespace LHCb { } // Compute chi2, note: we're using all hits (no crossed/uncrossed distinction) - std::tuple - MuonChamberMonitor::computeChi2( const GeomCache& geometryinfo, const xy_t& ext, const xy_t& extMS, - const CommonMuonHitRange& hitContainer ) const { + std::optional + LHCb::MuonChamberMonitor::computeChi2( const GeomCache& geometryinfo, const xy_t& ext, const xy_t& extMS2, + const CommonMuonHitRange& hitContainer ) const { - if ( hitContainer.size() == 0 ) return std::tuple( CommonMuonHit(), m_chi2bad, m_chi2bad ); // assign bad chi2 + if ( hitContainer.size() == 0 ) return chi2Result_t{ CommonMuonHit(), m_chi2bad, m_chi2bad }; // assign bad chi2 // take the hit with smallest chi2 assert( hitContainer.begin() != hitContainer.end() ); const auto chit = *std::min_element( hitContainer.begin(), hitContainer.end(), - [geometryinfo, ext, extMS]( const CommonMuonHit& hit1, const CommonMuonHit& hit2 ) { + [geometryinfo, ext, extMS2]( const CommonMuonHit& hit1, const CommonMuonHit& hit2 ) { auto dr2 = [&]( CommonMuonHit const& h ) { return pow( ( h.x() - ext.first ), 2 ) / - ( pow( geometryinfo.m_sigmapadX[h.station()][h.region()], 2 ) + extMS.first ) + + ( pow( geometryinfo.m_sigmapadX[h.station()][h.region()], 2 ) + extMS2.first ) + pow( ( h.y() - ext.second ), 2 ) / - ( pow( geometryinfo.m_sigmapadY[h.station()][h.region()], 2 ) + extMS.second ); + ( pow( geometryinfo.m_sigmapadY[h.station()][h.region()], 2 ) + extMS2.second ); }; return dr2( hit1 ) < dr2( hit2 ); } ); double chi2X = pow( ( chit.x() - ext.first ), 2 ) / - ( pow( geometryinfo.m_sigmapadX[chit.station()][chit.region()], 2 ) + extMS.first ); + ( pow( geometryinfo.m_sigmapadX[chit.station()][chit.region()], 2 ) + extMS2.first ); double chi2Y = pow( ( chit.y() - ext.second ), 2 ) / - ( pow( geometryinfo.m_sigmapadY[chit.station()][chit.region()], 2 ) + extMS.second ); - - return std::tuple( chit, chi2X, chi2Y ); + ( pow( geometryinfo.m_sigmapadY[chit.station()][chit.region()], 2 ) + extMS2.second ); + double winX = n_chi2cut * ( pow( geometryinfo.m_sigmapadX[chit.station()][chit.region()], 2 ) + extMS2.first ); + double winY = n_chi2cut * ( pow( geometryinfo.m_sigmapadY[chit.station()][chit.region()], 2 ) + extMS2.second ); + bool out_inner_acceptance = ( fabs( ext.first ) - sqrt( winX ) < geometryinfo.m_stationInnerX[chit.station()] ) && + ( fabs( ext.second ) - sqrt( winY ) < geometryinfo.m_stationInnerY[chit.station()] ); + bool out_outer_acceptance = ( fabs( ext.first ) + sqrt( winX ) > geometryinfo.m_stationOuterX[chit.station()] ) || + ( fabs( ext.second ) + sqrt( winY ) > geometryinfo.m_stationOuterY[chit.station()] ); + bool out_of_acceptance = out_inner_acceptance || out_outer_acceptance; + + return out_of_acceptance ? std::optional{} + : std::optional{ chi2Result_t{ chit, chi2X, chi2Y } }; } - int MuonChamberMonitor::HitsInChi2( const GeomCache& geometryinfo, const xy_t& ext, const xy_t& extMS, + int MuonChamberMonitor::HitsInChi2( const GeomCache& geometryinfo, const xy_t& ext, const xy_t& extMS2, const CommonMuonHitRange& hitContainer, int sign ) const { if ( hitContainer.size() == 0 ) return 0; int nhits = 0; for ( auto& hit : hitContainer ) { double chi2 = pow( ( hit.x() - sign * ext.first ), 2 ) / - ( pow( geometryinfo.m_sigmapadX[hit.station()][hit.region()], 2 ) + extMS.first ) + + ( pow( geometryinfo.m_sigmapadX[hit.station()][hit.region()], 2 ) + extMS2.first ) + pow( ( hit.y() - sign * ext.second ), 2 ) / - ( pow( geometryinfo.m_sigmapadY[hit.station()][hit.region()], 2 ) + extMS.second ); + ( pow( geometryinfo.m_sigmapadY[hit.station()][hit.region()], 2 ) + extMS2.second ); if ( chi2 < n_chi2cut ) { nhits++; } } return nhits;