diff --git a/Muon/MuonInterfaces/include/MuonInterfaces/MuonCluster.h b/Muon/MuonInterfaces/include/MuonInterfaces/MuonCluster.h
index 21df944cbab809837b4f3d7b894ac0a469f3fe4a..d25b817133117b2323d4921faf80d3484f97c838 100644
--- a/Muon/MuonInterfaces/include/MuonInterfaces/MuonCluster.h
+++ b/Muon/MuonInterfaces/include/MuonInterfaces/MuonCluster.h
@@ -94,21 +94,28 @@ public:
/// return the MuonTileIDs of the logical pads used for this cluster
std::vector<LHCb::Detector::Muon::TileID> getPadTiles() const;
/// number of logical pads in this hit
- int npads() const { return m_pads.size(); }
- int clsizeX() const { return m_xsize; } /// cluster size in X
- int clsizeY() const { return m_ysize; } /// cluster size in Y
- double dX() const { return m_dx; } /// error on X position
- double dY() const { return m_dy; } /// error on Y position
- double dZ() const { return m_dz; } /// error on Z position
- double minX() const { return m_hit_minx; }
- double maxX() const { return m_hit_maxx; }
- double minY() const { return m_hit_miny; }
- double maxY() const { return m_hit_maxy; }
- double minZ() const { return m_hit_minz; }
- double maxZ() const { return m_hit_maxz; }
- double X() const { return m_position.X(); }
- double Y() const { return m_position.Y(); }
- double Z() const { return m_position.Z(); }
+ int npads() const { return m_pads.size(); }
+ int clsizeX() const { return m_xsize; } /// cluster size in X
+ int clsizeY() const { return m_ysize; } /// cluster size in Y
+ double dx() const { return m_dx; } /// error on X position
+ double dX() const { return m_dx; } /// error on X position
+ double dy() const { return m_dy; } /// error on Y position
+ double dY() const { return m_dy; } /// error on Y position
+ double dz() const { return m_dz; } /// error on Z position
+ double dZ() const { return m_dz; } /// error on Z position
+ double minX() const { return m_hit_minx; }
+ double maxX() const { return m_hit_maxx; }
+ double minY() const { return m_hit_miny; }
+ double maxY() const { return m_hit_maxy; }
+ double minZ() const { return m_hit_minz; }
+ double maxZ() const { return m_hit_maxz; }
+ double x() const { return m_position.X(); }
+ double X() const { return m_position.X(); }
+ double y() const { return m_position.Y(); }
+ double Y() const { return m_position.Y(); }
+ double z() const { return m_position.Z(); }
+ double Z() const { return m_position.Z(); }
+
Gaudi::XYZVector pos() const { return m_position; }
MuonCluster& SetXYZ( double x, double y, double z ) {
m_position = { x, y, z };
diff --git a/Muon/MuonTools/CMakeLists.txt b/Muon/MuonTools/CMakeLists.txt
index 4b32be5ba62246359353efe9773f3b1ea2349161..856eec278b1283da4e4747f6fc71148ccc4a821c 100644
--- a/Muon/MuonTools/CMakeLists.txt
+++ b/Muon/MuonTools/CMakeLists.txt
@@ -21,6 +21,7 @@ gaudi_add_module(MuonTools
src/MuonPIDV2ToMuonTracks.cpp
src/MergeMuonPIDs.cpp
src/MergeMuonPIDsV1.cpp
+ src/MuonHitContainerToCommonMuonHits.cpp
LINK
Gaudi::GaudiAlgLib
Gaudi::GaudiKernel
diff --git a/Muon/MuonTools/src/MuonHitContainerToCommonMuonHits.cpp b/Muon/MuonTools/src/MuonHitContainerToCommonMuonHits.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..46eb1d5197e87cd283ce1386237a925303a9e4eb
--- /dev/null
+++ b/Muon/MuonTools/src/MuonHitContainerToCommonMuonHits.cpp
@@ -0,0 +1,52 @@
+/*****************************************************************************\
+* (c) Copyright 2000-2018 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. *
+\*****************************************************************************/
+
+// Gaudi
+#include "LHCbAlgs/Transformer.h"
+
+// LHCb
+#include "Event/PrHits.h"
+#include "MuonDAQ/CommonMuonHit.h"
+#include "MuonDAQ/MuonHitContainer.h"
+#include <Detector/Muon/MuonConstants.h>
+
+namespace LHCb {
+
+ /**
+ * Converter from MuonPIDs to MuonTracks so that alignment can work
+ *
+ * This has been developed as a work around after the replacement of MuonIDAlgLite by
+ * MuonIDHlt2Alg, as the former one was creating a list of MuonTracks while the later
+ * does not.
+ */
+
+ using CommonMuonHits = std::vector<CommonMuonHit, LHCb::Allocators::EventLocal<CommonMuonHit>>;
+
+ class MuonHitContainerToCommonMuonHits : public Algorithm::Transformer<CommonMuonHits( const MuonHitContainer& )> {
+
+ public:
+ MuonHitContainerToCommonMuonHits( const std::string& name, ISvcLocator* pSvcLocator )
+ : Transformer( name, pSvcLocator, KeyValue{ "Input", "" }, KeyValue{ "Output", "" } ) {}
+
+ CommonMuonHits operator()( const MuonHitContainer& hitContainer ) const override {
+ CommonMuonHits out;
+ out.reserve( hitContainer.hits( 0 ).size() * LHCb::Detector::Muon::nStations );
+ for ( std::size_t s = 0; s < LHCb::Detector::Muon::nStations; ++s ) {
+ auto const& hits = hitContainer.hits( s );
+ out.insert( out.end(), hits.begin(), hits.end() );
+ }
+ return out;
+ }
+ };
+
+ DECLARE_COMPONENT_WITH_ID( MuonHitContainerToCommonMuonHits, "MuonHitContainerToCommonMuonHits" )
+
+} // namespace LHCb
diff --git a/Tr/TrackTools/src/StandaloneMuonRec.cpp b/Tr/TrackTools/src/StandaloneMuonRec.cpp
index 7d40d979d5ecc5c6e18fcc1de82336342c209540..5372c690203e141c62cb2a049ab616dace84b5f6 100644
--- a/Tr/TrackTools/src/StandaloneMuonRec.cpp
+++ b/Tr/TrackTools/src/StandaloneMuonRec.cpp
@@ -28,7 +28,9 @@
#include "Magnet/DeMagnet.h"
#include "MuonDet/DeMuonDetector.h"
#include "MuonDet/MuonNamespace.h"
+#include "MuonInterfaces/MuonCluster.h"
#include "vdt/vdtMath.h"
+#include <Detector/Muon/MuonConstants.h>
#include "StandaloneMuonTrack.h"
@@ -67,14 +69,14 @@ namespace {
60., 120., 240., 480. } }; // M4
// -- Set tolerances for hit search in region
- constexpr std::array<float, 4> m_tolForRegion{ { 2.0, 4.0, 8.0, 10.0 } };
+ constexpr std::array<float, LHCb::Detector::Muon::nRegions> m_tolForRegion{ { 2.0, 4.0, 8.0, 10.0 } };
class Cache {
public:
- std::array<float, 4> stationZ{};
- Gaudi::XYZVector bdl;
- double zCenter{};
- int fieldPolarity{};
+ std::array<float, LHCb::Detector::Muon::nStations> stationZ{};
+ Gaudi::XYZVector bdl;
+ double zCenter{};
+ int fieldPolarity{};
Cache(){};
Cache( DeMuonDetector const& det ) {
for ( int s = 0; s != det.stations(); ++s ) { stationZ[s] = det.getStationZ( s ); }
@@ -82,41 +84,53 @@ namespace {
};
} // namespace
-class StandaloneMuonRec : public LHCb::Algorithm::Transformer<LHCb::Tracks( const MuonHitContainer&, const Cache& ),
- LHCb::Algorithm::Traits::usesConditions<Cache>> {
+/// MuonPositionConcept makes sure we are dealing either with CommonMuonHit or with MuonCluster
+template <typename MuonPosition>
+concept MuonPositionConcept =
+ requires { requires std::is_same_v<MuonPosition, CommonMuonHit> || std::is_same_v<MuonPosition, MuonCluster>; };
+
+template <MuonPositionConcept MuonPosition>
+class StandaloneMuonRec
+ : public LHCb::Algorithm::Transformer<
+ LHCb::Tracks( const std::vector<MuonPosition, LHCb::Allocators::EventLocal<MuonPosition>>&, const Cache& ),
+ LHCb::Algorithm::Traits::usesConditions<Cache>> {
public:
- using base_class_t = LHCb::Algorithm::Transformer<LHCb::Tracks( const MuonHitContainer&, const Cache& ),
+ using MuonPositionContainer = std::vector<MuonPosition, LHCb::Allocators::EventLocal<MuonPosition>>;
+ using PositionsInStations = std::array<std::vector<MuonPosition>, LHCb::Detector::Muon::nStations>;
+
+ using base_class_t = LHCb::Algorithm::Transformer<LHCb::Tracks( const MuonPositionContainer&, const Cache& ),
LHCb::Algorithm::Traits::usesConditions<Cache>>;
using base_class_t::addConditionDerivation;
/// Standard constructor
StandaloneMuonRec( const std::string& name, ISvcLocator* pSvcLocator )
- : base_class_t( name, pSvcLocator,
- { KeyValue{ "MuonHitsLocation", MuonHitContainerLocation::Default },
- KeyValue{ "ConditionsCache", "StandaloneMuonAlg-" + name + "-ConditionsCache" } },
- KeyValue{ "OutputMuonTracks", "Rec/Track/Muon" } ) {}
+ : base_class_t(
+ name, pSvcLocator,
+ { typename base_class_t::KeyValue{ "MuonHitsLocation", MuonHitContainerLocation::Default },
+ typename base_class_t::KeyValue{ "ConditionsCache", "StandaloneMuonAlg-" + name + "-ConditionsCache" } },
+ typename base_class_t::KeyValue{ "OutputMuonTracks", "Rec/Track/Muon" } ) {}
StatusCode initialize() override {
return base_class_t::initialize().andThen( [&] {
- this->addConditionDerivation( { DeMuonLocation::Default, LHCb::Det::Magnet::det_path },
- this->inputLocation<Cache>(),
- [&]( const DeMuonDetector& det, const DeMagnet& magnet ) {
- Cache cache{ det };
- const Gaudi::XYZPoint begin( 0., 0., 0. );
- const Gaudi::XYZPoint end( 0., 0., cache.stationZ[M2] );
- m_bIntegrator->calculateBdlAndCenter( magnet.fieldGrid(), begin, end, 0.0001, 0.,
- cache.zCenter, cache.bdl );
- cache.fieldPolarity = magnet.isDown() ? -1 : 1;
- debug() << "Integrated B field is " << cache.bdl.x() << " Tm"
- << " centered at Z=" << cache.zCenter / Gaudi::Units::m << " m"
- << " with polarity " << cache.fieldPolarity << endmsg;
- return cache;
- } );
+ this->template addConditionDerivation(
+ { DeMuonLocation::Default, LHCb::Det::Magnet::det_path }, this->template inputLocation<Cache>(),
+ [&]( const DeMuonDetector& det, const DeMagnet& magnet ) {
+ Cache cache{ det };
+ const Gaudi::XYZPoint begin( 0., 0., 0. );
+ const Gaudi::XYZPoint end( 0., 0., cache.stationZ[M2] );
+ m_bIntegrator->calculateBdlAndCenter( magnet.fieldGrid(), begin, end, 0.0001, 0., cache.zCenter,
+ cache.bdl );
+ cache.fieldPolarity = magnet.isDown() ? -1 : 1;
+ this->debug() << "Integrated B field is " << cache.bdl.x() << " Tm"
+ << " centered at Z=" << cache.zCenter / Gaudi::Units::m << " m"
+ << " with polarity " << cache.fieldPolarity << endmsg;
+ return cache;
+ } );
} );
}
- LHCb::Tracks operator()( const MuonHitContainer& hitContainer, const Cache& cache ) const override;
+ LHCb::Tracks operator()( const MuonPositionContainer& hitContainer, const Cache& cache ) const override;
ToolHandle<IBIntegrator> m_bIntegrator{ this, "BIntegrator", "BIntegrator" };
@@ -130,15 +144,18 @@ private:
Gaudi::Property<std::vector<float>> m_ParabolicCorrection{ this, "ParabolicCorrection", { 1.04, 0.14 } };
Gaudi::Property<std::vector<float>> m_resParams{ this, "m_resParams", { 0.015, 0.29 } };
Gaudi::Property<float> m_Constant{ this, "ConstantCorrection", 0., "In MeV" };
+ Gaudi::Property<float> m_yAt0Error{ this, "yAt0Error", std::numeric_limits<float>::max(),
+ "Error on y constraint to point to 0/0/0" };
- std::vector<StandaloneMuonTrack> muonSearch( const MuonHitContainer& hitContainer, const Cache& cache ) const;
+ std::vector<StandaloneMuonTrack<MuonPosition>> muonSearch( const PositionsInStations& hitContainer,
+ const Cache& cache ) const;
bool findCoincidence( const float x, const float y, const unsigned int station, const unsigned int regionBefore,
- const CommonMuonHitRange& hits, CommonMuonHit& hitcand ) const;
- void findmuonTrack( const MuonHitContainer& hitContainer, const Cache& cache,
- std::array<CommonMuonHit, 4>& bestCandidates, const int seed,
- std::vector<StandaloneMuonTrack>& muonTracks ) const;
- void recMomentum( StandaloneMuonTrack& muonTrack, const Cache& cache, LHCb::Track& track ) const;
- void detectClone( std::vector<StandaloneMuonTrack>& muonTracks, const Cache& cache ) const;
+ const std::vector<MuonPosition>& hits, MuonPosition& hitcand ) const;
+ void findmuonTrack( const PositionsInStations& hitContainer, const Cache& cache,
+ std::array<MuonPosition, LHCb::Detector::Muon::nStations>& bestCandidates, const int seed,
+ std::vector<StandaloneMuonTrack<MuonPosition>>& muonTracks ) const;
+ void recMomentum( StandaloneMuonTrack<MuonPosition>& muonTrack, const Cache& cache, LHCb::Track& track ) const;
+ void detectClone( std::vector<StandaloneMuonTrack<MuonPosition>>& muonTracks, const Cache& cache ) const;
// counters
mutable Gaudi::Accumulators::Counter<> m_countEvents{ this, "nEvents" };
@@ -148,26 +165,33 @@ private:
mutable Gaudi::Accumulators::MsgCounter<MSG::INFO> m_error_zeroBint{ this, "B integral is 0!!" };
};
-DECLARE_COMPONENT( StandaloneMuonRec )
+DECLARE_COMPONENT_WITH_ID( StandaloneMuonRec<CommonMuonHit>, "StandaloneMuonRec" )
+DECLARE_COMPONENT_WITH_ID( StandaloneMuonRec<MuonCluster>, "StandaloneMuonRecWithClusters" )
//=============================================================================
// Main execution
//=============================================================================
-LHCb::Tracks StandaloneMuonRec::operator()( const MuonHitContainer& hitContainer, const Cache& cache ) const {
+template <MuonPositionConcept MuonPosition>
+LHCb::Tracks StandaloneMuonRec<MuonPosition>::operator()( const MuonPositionContainer& hitContainer,
+ const Cache& cache ) const {
- if ( msgLevel( MSG::DEBUG ) ) debug() << "==> Execute" << endmsg;
+ if ( this->msgLevel( MSG::DEBUG ) ) this->debug() << "==> Execute" << endmsg;
+
+ PositionsInStations posInStations;
+ for ( auto hit : hitContainer ) { posInStations[hit.station()].push_back( hit ); }
LHCb::Tracks outputTracks;
- outputTracks.reserve( 1024 );
+ outputTracks.reserve( 10 );
++m_countEvents;
- auto muonTracks = muonSearch( hitContainer, cache );
+ auto muonTracks = muonSearch( posInStations, cache );
if ( m_cloneKiller ) { detectClone( muonTracks, cache ); }
for ( auto& muTrack : muonTracks ) {
if ( muTrack.isClone() ) continue;
+ muTrack.setYAt0Err( m_yAt0Error );
auto sc = muTrack.linearFit();
if ( !sc ) {
++m_failed_linearfit;
@@ -181,14 +205,15 @@ LHCb::Tracks StandaloneMuonRec::operator()( const MuonHitContainer& hitContainer
}
m_countMuCandidates += outputTracks.size();
- if ( msgLevel( MSG::DEBUG ) ) debug() << " stored candidates " << outputTracks.size() << endmsg;
+ if ( this->msgLevel( MSG::DEBUG ) ) this->debug() << " stored candidates " << outputTracks.size() << endmsg;
return outputTracks;
}
-
-bool StandaloneMuonRec::findCoincidence( const float x, const float y, const unsigned int station,
- const unsigned int regionBefore, const CommonMuonHitRange& hits,
- CommonMuonHit& hitcand ) const {
+template <MuonPositionConcept MuonPosition>
+bool StandaloneMuonRec<MuonPosition>::findCoincidence( const float x, const float y, const unsigned int station,
+ const unsigned int regionBefore,
+ const std::vector<MuonPosition>& hits,
+ MuonPosition& hitcand ) const {
const auto tol = m_tolForRegion[regionBefore];
float deltaYmin = 9999.;
@@ -198,7 +223,8 @@ bool StandaloneMuonRec::findCoincidence( const float x, const float y, const uns
float deltaX = fabs( x - hit.x() );
float deltaY = fabs( y - hit.y() );
//-- Check if the hit is within the FOI (copy from MuonCombRec); in any case a xFOI >1000 is not considered
- if ( deltaX < m_Xmax[station * 4 + regionBefore] && deltaY < m_Ymax[station * 4 + regionBefore] &&
+ if ( deltaX < m_Xmax[station * LHCb::Detector::Muon::nStations + regionBefore] &&
+ deltaY < m_Ymax[station * LHCb::Detector::Muon::nStations + regionBefore] &&
( deltaY < deltaYmin - tol ||
( deltaY < deltaYmin + tol && ( deltaX < deltaXmin - tol || fabs( deltaXmin - deltaX ) < 0.1 ) ) ) ) {
deltaXmin = deltaX;
@@ -209,9 +235,11 @@ bool StandaloneMuonRec::findCoincidence( const float x, const float y, const uns
}
return findCand;
}
-void StandaloneMuonRec::findmuonTrack( const MuonHitContainer& hitContainer, const Cache& cache,
- std::array<CommonMuonHit, 4>& bestCandidates, const int seed,
- std::vector<StandaloneMuonTrack>& muonTracks ) const {
+template <MuonPositionConcept MuonPosition>
+void StandaloneMuonRec<MuonPosition>::findmuonTrack(
+ const PositionsInStations& hitContainer, const Cache& cache,
+ std::array<MuonPosition, LHCb::Detector::Muon::nStations>& bestCandidates, const int seed,
+ std::vector<StandaloneMuonTrack<MuonPosition>>& muonTracks ) const {
float xseed = bestCandidates[seed].x();
float yseed = bestCandidates[seed].y();
unsigned int hit_region = bestCandidates[seed].region();
@@ -220,7 +248,7 @@ void StandaloneMuonRec::findmuonTrack( const MuonHitContainer& hitContainer, con
bool findAll = false;
for ( int ista = seed - 1; ista > -1; ista-- ) {
- auto findCandidate = findCoincidence( x, y, ista, hit_region, hitContainer.hits( ista ), bestCandidates[ista] );
+ auto findCandidate = findCoincidence( x, y, ista, hit_region, hitContainer.at( ista ), bestCandidates[ista] );
if ( !findCandidate && ( !m_secondLoop || ista < M4 ) ) break;
if ( findCandidate && ista == M2 ) {
findAll = true;
@@ -243,34 +271,33 @@ void StandaloneMuonRec::findmuonTrack( const MuonHitContainer& hitContainer, con
}
if ( findAll ) {
// create the muon track
- StandaloneMuonTrack muon;
- muon.setPoint( 0, bestCandidates[M2] );
- muon.setPoint( 1, bestCandidates[M3] );
+ StandaloneMuonTrack<MuonPosition> muon;
+ muon.addPoint( bestCandidates[M2] );
+ muon.addPoint( bestCandidates[M3] );
if ( m_secondLoop && bestCandidates[M4].station() == bestCandidates[M5].station() ) {
- muon.setPoint( 2, bestCandidates[M4] );
- muon.setnHits( 3 );
+ muon.addPoint( bestCandidates[M4] );
} else {
- muon.setPoint( 2, bestCandidates[M4] );
- muon.setPoint( 3, bestCandidates[M5] );
- muon.setnHits( 4 );
+ muon.addPoint( bestCandidates[M4] );
+ muon.addPoint( bestCandidates[M5] );
}
muonTracks.push_back( muon );
}
}
-std::vector<StandaloneMuonTrack> StandaloneMuonRec::muonSearch( const MuonHitContainer& hitContainer,
- const Cache& cache ) const {
- std::vector<StandaloneMuonTrack> muonTracks;
- muonTracks.reserve( 48 );
-
- const auto& hitsM5 = hitContainer.hits( M5 );
- std::array<CommonMuonHit, 4> bestCandidates;
+template <MuonPositionConcept MuonPosition>
+std::vector<StandaloneMuonTrack<MuonPosition>>
+StandaloneMuonRec<MuonPosition>::muonSearch( const PositionsInStations& hitContainer, const Cache& cache ) const {
+ std::vector<StandaloneMuonTrack<MuonPosition>> muonTracks;
+ muonTracks.reserve( 12 );
+
+ const auto& hitsM5 = hitContainer.at( M5 );
+ std::array<MuonPosition, LHCb::Detector::Muon::nStations> bestCandidates;
for ( auto& hit : hitsM5 ) {
bestCandidates[3] = hit;
findmuonTrack( hitContainer, cache, bestCandidates, M5, muonTracks );
}
///---second loop from M4
if ( m_secondLoop ) {
- const auto& hitsM4 = hitContainer.hits( M4 );
+ const auto& hitsM4 = hitContainer[M4];
for ( auto& hit : hitsM4 ) {
// To remove these hits of M4 used in the first round of search
auto used = std::find_if( muonTracks.begin(), muonTracks.end(),
@@ -285,7 +312,9 @@ std::vector<StandaloneMuonTrack> StandaloneMuonRec::muonSearch( const MuonHitCon
}
// estimate the momentum of muonTrack
-void StandaloneMuonRec::recMomentum( StandaloneMuonTrack& track, const Cache& cache, LHCb::Track& lbtrack ) const {
+template <MuonPositionConcept MuonPosition>
+void StandaloneMuonRec<MuonPosition>::recMomentum( StandaloneMuonTrack<MuonPosition>& track, const Cache& cache,
+ LHCb::Track& lbtrack ) const {
const float bdlX = cache.bdl.x();
const float fieldPolarity = cache.fieldPolarity;
@@ -296,7 +325,7 @@ void StandaloneMuonRec::recMomentum( StandaloneMuonTrack& track, const Cache& ca
LHCb::State state( LHCb::StateVector( trackPos, Gaudi::XYZVector( track.sx(), track.sy(), 1.0 ), 1. / 10000. ) );
const auto Zend = cache.stationZ[M5];
- Gaudi::XYZPoint endtrackPos( track.bx() + track.sx() * Zfirst, track.by() + track.sy() * Zend, Zend );
+ Gaudi::XYZPoint endtrackPos( track.bx() + track.sx() * Zend, track.by() + track.sy() * Zend, Zend );
LHCb::State endstate(
LHCb::StateVector( endtrackPos, Gaudi::XYZVector( track.sx(), track.sy(), 1.0 ), 1. / 10000. ) );
@@ -353,25 +382,38 @@ void StandaloneMuonRec::recMomentum( StandaloneMuonTrack& track, const Cache& ca
state.setLocation( LHCb::State::Location::Muon );
endstate.setLocation( LHCb::State::Location::LastMeasurement );
- debug() << "Muon state = " << state << endmsg;
+ this->debug() << "Muon state = " << state << endmsg;
lbtrack.clearStates();
lbtrack.addToStates( state );
lbtrack.addToStates( endstate );
lbtrack.setChi2PerDoF( track.chi2x() + track.chi2y() );
- lbtrack.setNDoF( track.nHits() - 2 );
+ const int DoFForPointAtZero = m_yAt0Error < std::numeric_limits<float>::max();
+ lbtrack.setNDoF( track.nHits() - 2 + DoFForPointAtZero ); // add one degree of freedom for the point at 0/0/0
for ( int i = 0; i < track.nHits(); i++ ) {
- const auto Tile = track.point( i ).tile();
- lbtrack.addToLhcbIDs( ( LHCb::LHCbID )( Tile ) );
- debug() << " Muon Hit " << i << " tile " << Tile << " tiles in station " << track.point( i ).station() << endmsg;
+ if constexpr ( std::is_same_v<MuonPosition, CommonMuonHit> ) {
+ const auto tile = track.point( i ).tile();
+ lbtrack.addToLhcbIDs( ( LHCb::LHCbID )( tile ) );
+ this->debug() << " Muon Hit " << i << " tile " << tile << " tiles in station " << track.point( i ).station()
+ << endmsg;
+ } else {
+ const auto tiles = track.point( i ).getPadTiles();
+ for ( const auto tile : tiles ) {
+ lbtrack.addToLhcbIDs( ( LHCb::LHCbID )( tile ) );
+ this->debug() << " Muon Hit " << i << " tile " << tile << " tiles in station " << track.point( i ).station()
+ << endmsg;
+ }
+ }
}
lbtrack.setPatRecStatus( LHCb::Track::PatRecStatus::PatRecIDs );
lbtrack.setType( LHCb::Track::Types::Muon );
}
-void StandaloneMuonRec::detectClone( std::vector<StandaloneMuonTrack>& muonTracks, const Cache& cache ) const {
+template <MuonPositionConcept MuonPosition>
+void StandaloneMuonRec<MuonPosition>::detectClone( std::vector<StandaloneMuonTrack<MuonPosition>>& muonTracks,
+ const Cache& cache ) const {
for ( auto itMuonTrackFirst = muonTracks.begin(); itMuonTrackFirst < muonTracks.end(); itMuonTrackFirst++ ) {
for ( auto itMuonTrackSecond = itMuonTrackFirst + 1; itMuonTrackSecond < muonTracks.end(); itMuonTrackSecond++ ) {
diff --git a/Tr/TrackTools/src/StandaloneMuonTrack.h b/Tr/TrackTools/src/StandaloneMuonTrack.h
index 77df068d551b160df4b2b91d424ae1084f29a521..1e6c92c531a68bfc4bf48a0fad40e652c3249d88 100644
--- a/Tr/TrackTools/src/StandaloneMuonTrack.h
+++ b/Tr/TrackTools/src/StandaloneMuonTrack.h
@@ -13,9 +13,11 @@
#include "Core/FloatComparison.h"
#include "MuonDAQ/CommonMuonHit.h"
+#include "Detector/Muon/MuonConstants.h"
#include "GaudiKernel/Point3DTypes.h"
#include "GaudiKernel/Vector3DTypes.h"
+#include <boost/container/small_vector.hpp>
#include <string>
/** @class StandaloneMuonTrack StandaloneMuonTrack.h
@@ -29,18 +31,33 @@
* @date 2011-03-03
*/
namespace {
- void LinearFitXZ( bool XZ, std::array<CommonMuonHit, 4>& m_points, int nHits, float& slope, float& trunc,
- float& chi2ndof, float& err_slope, float& err_trunc, float& cov ) {
+
+ enum class FitMode { XZ, YZ };
+
+ struct FitResult {
+ float slope = -1e6;
+ float trunc = -1e6;
+ float chi2ndof = -1e6;
+ float err_slope = -1e6;
+ float err_trunc = -1e6;
+ float cov = -1e6;
+ };
+
+ template <typename MuonPositionClass, FitMode fitMode>
+ FitResult linearFitImpl( boost::container::small_vector<MuonPositionClass, LHCb::Detector::Muon::nStations>& points,
+ float yAt0Err = std::numeric_limits<float>::max() ) {
float sz2, sz, s0, sxz, sx, sx2;
sz2 = sz = s0 = sxz = sx = sx2 = 0;
- for ( int i = 0; i < nHits; i++ ) {
- const auto hit = m_points[i];
- const auto z = hit.z();
- auto x = hit.x();
- auto xerr = 2. * hit.dx();
- if ( !XZ ) {
- x = hit.y();
- xerr = 2 * hit.dy();
+
+ FitResult result;
+
+ for ( const auto& point : points ) {
+ const auto z = point.z();
+ auto x = point.x();
+ auto xerr = 2. * point.dx();
+ if ( fitMode == FitMode::YZ ) {
+ x = point.y();
+ xerr = 2 * point.dy();
}
sz2 += z * z / xerr / xerr;
sz += z / xerr / xerr;
@@ -49,91 +66,86 @@ namespace {
sx += x / xerr / xerr;
sx2 += x * x / xerr / xerr;
}
- float xdet = sz2 * s0 - sz * sz;
- if ( !LHCb::essentiallyZero( xdet ) ) {
- slope = ( sxz * s0 - sx * sz ) / xdet;
- trunc = ( sx * sz2 - sxz * sz ) / xdet;
- err_trunc = sqrt( sz2 / xdet );
- err_slope = sqrt( s0 / xdet );
- cov = -sz / xdet;
+ if ( yAt0Err < std::numeric_limits<float>::max() )
+ s0 += 1.0 / ( yAt0Err * yAt0Err ); // this applies a constraint to 0/0/0
- chi2ndof = ( sx2 + slope * slope * sz2 + trunc * trunc * s0 - 2. * slope * sxz - 2. * trunc * sx +
- 2 * slope * trunc * sz ) /
- ( nHits - 2 );
+ const float xdet = sz2 * s0 - sz * sz;
+ if ( !LHCb::essentiallyZero( xdet ) ) {
+ result.slope = ( sxz * s0 - sx * sz ) / xdet;
+ result.trunc = ( sx * sz2 - sxz * sz ) / xdet;
+
+ result.err_trunc = sqrt( sz2 / xdet );
+ result.err_slope = sqrt( s0 / xdet );
+ result.cov = -sz / xdet;
+
+ result.chi2ndof =
+ ( sx2 + result.slope * result.slope * sz2 + result.trunc * result.trunc * s0 - 2. * result.slope * sxz -
+ 2. * result.trunc * sx + 2 * result.slope * result.trunc * sz ) /
+ ( points.size() - 2 + ( yAt0Err < std::numeric_limits<float>::max() ) ); // add one degree of freedom if we
+ // add the "point" at 0/0/0
}
+
+ return result;
}
} // namespace
+template <typename MuonPositionClass>
class StandaloneMuonTrack final {
public:
/// Standard constructor
- StandaloneMuonTrack() { m_clone = 0; };
+ StandaloneMuonTrack() {
+ m_clone = false;
+ m_points.reserve( 4 );
+ };
virtual ~StandaloneMuonTrack(){}; ///< Destructor
- void setPoint( unsigned int station, CommonMuonHit point ) { m_points[station] = point; };
+ void setPoint( unsigned int station, MuonPositionClass point ) { m_points[station] = point; };
+ void addPoint( MuonPositionClass point ) { m_points.push_back( point ); };
- CommonMuonHit point( unsigned int station ) const { return m_points[station]; };
+ MuonPositionClass point( unsigned int station ) const { return m_points[station]; };
- void setClone() { m_clone = 1; };
- bool isClone() const { return m_clone > 0; }
+ void setClone() { m_clone = true; };
+ bool isClone() const { return m_clone; }
- double slopeX( int stationFirst, int stationSecond, double zFirst, double zSecond ) const {
- return ( m_points[stationFirst].x() - m_points[stationSecond].x() ) / ( zFirst - zSecond );
- };
- double slopeY( int stationFirst, int stationSecond, double zFirst, double zSecond ) const {
- return ( m_points[stationFirst].y() - m_points[stationSecond].y() ) / ( zFirst - zSecond );
- };
+ // -- used to constrain the straight line fit in y to 0/0/0
+ void setYAt0Err( const float yAt0Err ) { m_yAt0Err = yAt0Err; }
// Fit with a min chi^2 in the 2 projections xz and yz
bool linearFit() {
- double dof = nHits() - 2.;
+ const double dof = nHits() - 2.;
if ( dof < 0 ) return false;
- LinearFitXZ( true, m_points, nHits(), m_sx, m_bx, m_chi2x, m_errsx, m_errbx, m_covbsx );
- LinearFitXZ( false, m_points, nHits(), m_sy, m_by, m_chi2y, m_errsy, m_errby, m_covbsy );
- if ( m_chi2x > -1.f && m_chi2y > -1.f )
- return true;
- else
- return false;
+ m_resultXZ = linearFitImpl<MuonPositionClass, FitMode::XZ>( m_points );
+ m_resultYZ = linearFitImpl<MuonPositionClass, FitMode::YZ>( m_points, m_yAt0Err );
+ return m_resultXZ.chi2ndof > -1.f && m_resultYZ.chi2ndof > -1.f;
};
- inline double chi2x() const { return m_chi2x; } /// chi2/dof XZ
- inline double chi2y() const { return m_chi2y; } /// chi2/dof YZ
+ inline double chi2x() const { return m_resultXZ.chi2ndof; } /// chi2/dof XZ
+ inline double chi2y() const { return m_resultYZ.chi2ndof; } /// chi2/dof YZ
// slope XZ
- inline double sx() const { return m_sx; }
+ inline double sx() const { return m_resultXZ.slope; }
// intercept XZ
- inline double bx() const { return m_bx; }
+ inline double bx() const { return m_resultXZ.trunc; }
// slope YZ
- inline double sy() const { return m_sy; }
+ inline double sy() const { return m_resultYZ.slope; }
// intercept YZ
- inline double by() const { return m_by; }
+ inline double by() const { return m_resultYZ.trunc; }
// errors on the above parameters
- inline double errsx() const { return m_errsx; }
- inline double errbx() const { return m_errbx; }
- inline double covbsx() const { return m_covbsx; }
- inline double errsy() const { return m_errsy; }
- inline double errby() const { return m_errby; }
- inline double covbsy() const { return m_covbsy; }
+ inline double errsx() const { return m_resultXZ.err_slope; }
+ inline double errbx() const { return m_resultXZ.err_trunc; }
+ inline double covbsx() const { return m_resultXZ.cov; }
+ inline double errsy() const { return m_resultYZ.err_slope; }
+ inline double errby() const { return m_resultYZ.err_trunc; }
+ inline double covbsy() const { return m_resultYZ.cov; }
- inline void setnHits( int n ) { m_nHits = n; }
- inline int nHits() const { return m_nHits; }
+ inline int nHits() const { return m_points.size(); }
private:
- std::array<CommonMuonHit, 4> m_points;
- unsigned int m_clone;
- unsigned int m_nHits;
-
- float m_chi2x = -1.f;
- float m_chi2y = -1.f;
- float m_sx = -1.f;
- float m_bx = -1.f;
- float m_sy = -1.f;
- float m_by = -1.f;
- float m_errsx = -1.f;
- float m_errbx = -1.f;
- float m_covbsx = -1.f;
- float m_errsy = -1.f;
- float m_errby = -1.f;
- float m_covbsy = -1.f;
+ boost::container::small_vector<MuonPositionClass, LHCb::Detector::Muon::nStations> m_points;
+ bool m_clone;
+ float m_yAt0Err = std::numeric_limits<float>::max();
+
+ FitResult m_resultXZ;
+ FitResult m_resultYZ;
};