From e539800a85beabe41885d40140bce45d0779d3d4 Mon Sep 17 00:00:00 2001 From: Olli Lupton <oliver.lupton@cern.ch> Date: Tue, 14 Nov 2017 16:42:52 +0100 Subject: [PATCH 1/5] Various tweaks to get a standalone muon reconstruction running in the upgrade (where M1 is removed, and .station() == 0 means M2...) --- .../src/component/MuonCombRec.cpp | 64 ++-- Muon/MuonTrackRec/src/component/MuonCombRec.h | 39 +-- .../src/component/MuonPadFromCoord.cpp | 46 ++- .../src/component/MuonPadFromCoord.h | 12 +- .../src/component/MuonTrackMomRec.cpp | 5 +- Tr/TrackTools/src/MuonSeeding.cpp | 291 ++++++++++++++++++ Tr/TrackTools/src/MuonSeeding.h | 76 +++++ 7 files changed, 468 insertions(+), 65 deletions(-) create mode 100644 Tr/TrackTools/src/MuonSeeding.cpp create mode 100644 Tr/TrackTools/src/MuonSeeding.h diff --git a/Muon/MuonTrackRec/src/component/MuonCombRec.cpp b/Muon/MuonTrackRec/src/component/MuonCombRec.cpp index 59d5821f3f7..2ff18db7bc6 100644 --- a/Muon/MuonTrackRec/src/component/MuonCombRec.cpp +++ b/Muon/MuonTrackRec/src/component/MuonCombRec.cpp @@ -38,7 +38,7 @@ DECLARE_COMPONENT( MuonCombRec ) namespace { // -- initialize the pad size. Hardwired to speed up. - constexpr std::array<double, 20> m_padX = {{ + constexpr std::array<double, 20> m_padX = {{ // R1 R2 R3 R4 10.1667 , 20.1667 , 40.333 , 80.6667 , // M1 6.41667, 12.75 , 25.5 , 51. , // M2 @@ -47,7 +47,7 @@ namespace { 31.6667 , 63. , 126. , 252 }}; // M5 // R1 R2 R3 R4 - constexpr std::array<double, 20> m_padY = {{ + constexpr std::array<double, 20> m_padY = {{ // R1 R2 R3 R4 24.7586 , 49.7584 , 99.7581 , 199.757 , // M1 31.3585 , 62.9583 , 126.158 , 252.557 , // M2 @@ -157,9 +157,18 @@ StatusCode MuonCombRec::initialize() { // incident service incSvc()->addListener( this, IncidentType::EndEvent ); + // OL 171114: offset introduced + // the idea is that the various arrays will always be indexed following the M1, M2 etc. enums + // geometry stuff ---> m_nStation = m_muonDetector->stations(); - if(m_nStation != 5){ + if ( m_nStation == 5 ) { + // Run 1 or 2 + m_stationOffset = 0; + } else if ( m_nStation == 4 ) { + // Upgrade + m_stationOffset = 1; + } else { err()<<"Wrong station number: "<<m_nStation<<endmsg; return StatusCode::FAILURE; } @@ -173,9 +182,9 @@ StatusCode MuonCombRec::initialize() { } // -- get the z position of stations - for (int station=0; station<m_nStation; ++station){ - m_zStations[station] = m_muonDetector->getStationZ(station); - if ( msgLevel(MSG::DEBUG) ) debug() << "Z of station M" << station+1 << ": " << m_zStations[station] << endmsg; + for (int station = 0; station < m_nStation; ++station){ + m_zStations[station + m_stationOffset] = m_muonDetector->getStationZ(station); + if ( msgLevel(MSG::DEBUG) ) debug() << "Z of station M" << station+m_stationOffset+1 << ": " << m_zStations[station+m_stationOffset] << endmsg; } // <--- @@ -196,15 +205,18 @@ StatusCode MuonCombRec::initialize() { // local stuff if ( msgLevel(MSG::DEBUG) ){ debug()<<"xFOIs: "<<endmsg; - for(int is = 0; is < m_nStation-1; is++){ + for(int is = m_stationOffset; is < m_nStation + m_stationOffset - 1; is++) { // [M1] M2-4 + debug() << "M" << is+1 << " "; for(int ir = 0; ir < m_nRegion; ir++){ int key = is*4+ir; debug()<< m_xFOIs[key] <<" "; } debug()<<endmsg; } + debug()<<"yFOIs: "<<endmsg; - for(int is = 0; is < m_nStation-1; is++){ + for(int is = m_stationOffset; is < m_nStation + m_stationOffset - 1; is++) { + debug() << "M" << is+1 << " "; for(int ir = 0; ir < m_nRegion; ir++){ int key = is*4+ir; debug()<< m_yFOIs[key] <<" "; @@ -215,7 +227,8 @@ StatusCode MuonCombRec::initialize() { if ( msgLevel(MSG::DEBUG) ){ debug()<< "x pad sizes: " << endmsg; - for(int is = 0; is < m_nStation; is++){ + for(int is = m_stationOffset; is < m_nStation + m_stationOffset; is++){ + debug() << "M" << is+1 << " "; for(int ir = 0; ir < m_nRegion; ir++){ int key = is*4+ir; debug()<< m_padX[key] << " "; @@ -223,7 +236,8 @@ StatusCode MuonCombRec::initialize() { debug() << endmsg; } debug()<< "y pad sizes: " << endmsg; - for(int is = 0; is < m_nStation; is++){ + for(int is = m_stationOffset; is < m_nStation + m_stationOffset; is++){ + debug() << "M" << is+1 << " "; for(int ir = 0; ir < m_nRegion; ir++){ int key = is*4+ir; debug()<< m_padY[key] << " "; @@ -389,9 +403,9 @@ StatusCode MuonCombRec::muonSearch() { // -- Local array x-talk hits std::vector<const MuonHit*> xtHits; - int SS = m_seedStation; // SS = Seed Station + int SS = m_seedStation; // SS = Seed Station (numbering scheme where 0 = M1, 1 = M2 etc. - if ( msgLevel(MSG::DEBUG) ) debug() << "Seed Station is: " << SS << endmsg; + if ( msgLevel(MSG::DEBUG) ) debug() << "Seed Station is: M" << SS+1 << endmsg; // -- Loop over SS hits and copy them in a local array with SS seeds for ( int region = 0; region < m_nRegion; ++region){ @@ -408,7 +422,7 @@ StatusCode MuonCombRec::muonSearch() { candidates[SS].clear(); candidates[SS].push_back(seedHit); - if ( msgLevel(MSG::DEBUG) ) debug() << "Seed : " << seedHit->station() << " " + if ( msgLevel(MSG::DEBUG) ) debug() << "Seed : M" << seedHit->station() + m_stationOffset + 1 << " " << seedHit->region() << " " << seedHit->X() << " " << seedHit->Y() << " " @@ -429,7 +443,7 @@ StatusCode MuonCombRec::muonSearch() { }); if(used!=xtHits.end()) { if ( msgLevel(MSG::DEBUG) ) debug() << " Hit already used for another seed " - << (*used)->station() << " " + << "M" << ((*used)->station() + m_stationOffset + 1) << " " << (*used)->region() << " " << (*used)->X() << " " << (*used)->Y() << " " @@ -451,11 +465,11 @@ StatusCode MuonCombRec::muonSearch() { const MuonHit* seedHit2 = *itSSn; - int istam = seedHit->station(); + int istam = seedHit->station() + m_stationOffset; // in upgrade station() returns 0 for M2 int iregm = seedHit->region(); int keym = istam*4+iregm; - int istan = seedHit2->station(); + int istan = seedHit2->station() + m_stationOffset; int iregn = seedHit2->region(); int keyn = istan*4+iregn; @@ -481,7 +495,7 @@ StatusCode MuonCombRec::muonSearch() { debug() << "Accumulated x-talk points for this seed: -->" << endmsg; for( const MuonHit* hit : seedXTHits ){ - debug() << hit->station() << " " << hit->region() + debug() << "M" << hit->station() + m_stationOffset + 1 << " " << hit->region() << " " << hit->X() << " " << hit->Y() << " " << hit->Z() <<endmsg; @@ -549,7 +563,7 @@ StatusCode MuonCombRec::muonSearch() { bool matchingFound = true; - for(int is = nextStation; is >= 0; --is){ + for(int is = nextStation; is >= m_stationOffset; --is){ if(is == m_skipStation) continue; // skip a station @@ -646,7 +660,7 @@ StatusCode MuonCombRec::muonSearch() { // match hits with the track extrapolation //======================================================================== StatusCode MuonCombRec::findMatching( const double x, const double y, - const int station, const int regionBefore, + const int station, const int regionBefore, // station follows the M1=0, M2=1 scheme std::vector<const MuonHit*> &candidates ){ double toll = m_tolForRegion[regionBefore]; @@ -716,11 +730,11 @@ StatusCode MuonCombRec::findMatching( const double x, const double y, if(deltaX<2&&deltaY<2) continue; // same pad...;-) - int csta = candidate->station(); + int csta = candidate->station() + m_stationOffset; int creg = candidate->region(); int ckey = csta*4+creg; - int ista = hit->station(); + int ista = hit->station() + m_stationOffset; int ireg = hit->region(); int ikey = ista*4+ireg; @@ -748,7 +762,7 @@ StatusCode MuonCombRec::cloneKiller() for(const auto& hitT1 : hitsT1) { // loop over hits of 1st trk - int st1 = hitT1->station(); + int st1 = hitT1->station() + m_stationOffset; if(st1 != M2 && st1 != M3 && st1 != M4) continue;//look only in M2,M3,M4 for(const auto& hitT2 : hitsT2){ // loop over hits of 2nd trk @@ -806,7 +820,7 @@ StatusCode MuonCombRec::strongCloneKiller() for(const auto& hitT1 : hitsT1){ // loop over hits of 1st trk - int st1 = hitT1->station(); + int st1 = hitT1->station() + m_stationOffset; if(st1 != M2 && st1 != M3) continue; // look only in M2,M3 for(const auto& hitT2 : hitsT2){ // loop over hits of 2nd trk @@ -855,7 +869,7 @@ StatusCode MuonCombRec::sortMuonHits(){ // -- Fill the sorted hit array for later track reconstruction for( const auto& hit : *m_trackhits ){ - int station = hit.station(); + int station = hit.station() + m_stationOffset; int region = hit.region(); m_sortedHits[station*4+region].push_back(&hit); } @@ -863,7 +877,7 @@ StatusCode MuonCombRec::sortMuonHits(){ if ( msgLevel(MSG::DEBUG) ){ for( const auto& hitVec : m_sortedHits){ for( const MuonHit* hit : hitVec ){ - debug() << "sorted hit stored station: " << hit->station() + debug() << "sorted hit stored station: M" << hit->station() + m_stationOffset + 1 <<" sorted hit stored region: "<< hit->region() <<" x,y,z: " << hit->X() << " " diff --git a/Muon/MuonTrackRec/src/component/MuonCombRec.h b/Muon/MuonTrackRec/src/component/MuonCombRec.h index 5600f80aa70..0d44f766c4d 100644 --- a/Muon/MuonTrackRec/src/component/MuonCombRec.h +++ b/Muon/MuonTrackRec/src/component/MuonCombRec.h @@ -141,7 +141,8 @@ private: void clear(); - int m_nStation; + int m_stationOffset; // offset into lookup tables (0 in Run1/2, 1 in upgrade) + int m_nStation; // number of stations (5 in Run1/2, 4 in upgrade) int m_nRegion; // -- timing @@ -156,48 +157,48 @@ private: // enable timers Gaudi::Property<bool> m_measureTime {this, "MeasureTime", false}; - + // enable clone finding and killing bool m_cloneKiller; - Gaudi::Property<bool> m_optCloneKiller - {this, "CloneKiller", true, + Gaudi::Property<bool> m_optCloneKiller + {this, "CloneKiller", true, "enable clone finding and killing"}; // enable strong clone finding and killing bool m_strongCloneKiller; - Gaudi::Property<bool> m_optStrongCloneKiller + Gaudi::Property<bool> m_optStrongCloneKiller {this, "StrongCloneKiller", true, "enable strong clone finding and killing"}; // station to be skipped (e.g. for eff. study): 0-4; // -1 = keep all stns (default) int m_skipStation; - Gaudi::Property<int> m_optSkipStation + Gaudi::Property<int> m_optSkipStation {this, "SkipStation", -1, "station to be skipped (e.g. for eff. study): 0-4; -1 = keep all stns (default)"}; - Gaudi::Property<bool> m_physicsTiming - {this, "PhysicsTiming", true, + Gaudi::Property<bool> m_physicsTiming + {this, "PhysicsTiming", true, "if true, we assume that timing is the final one (accounting for TOF from primary vx)"}; - Gaudi::Property<bool> m_assumeCosmics - {this, "AssumeCosmics", false, + Gaudi::Property<bool> m_assumeCosmics + {this, "AssumeCosmics", false, "if true (default) we assume that tracks are of cosmic origin (can be backward)"}; - Gaudi::Property<bool> m_assumePhysics - {this, "AssumePhysics", true, + Gaudi::Property<bool> m_assumePhysics + {this, "AssumePhysics", true, "if true we assume that tracks have the 'right' direction (pz>0)"}; Gaudi::Property<std::string> m_decToolName {this, "DecodingTool", "MuonHitDecode" "name of decoding tool (MuonHitDecode for offline, MuonMonHitDecode for online monitoring)"}; - Gaudi::Property<std::string> m_padToolName - {this, "PadRecTool", "MuonPadRec", + Gaudi::Property<std::string> m_padToolName + {this, "PadRecTool", "MuonPadRec", "name of pad rec tool (MuonPadRec only option so far)"}; - - Gaudi::Property<std::string> m_clusterToolName - {this, "ClusterTool", "MuonFakeClustering", + + Gaudi::Property<std::string> m_clusterToolName + {this, "ClusterTool", "MuonFakeClustering", "name of clustering tool"}; /// cross talk @@ -207,8 +208,8 @@ private: int m_seedStation; Gaudi::Property<int> m_optSeedStation {this, "SeedStation", 4, "default seed: M5"}; - Gaudi::Property<std::string> m_trackOutputLoc - {this, "TracksOutputLocation", "Rec/Muon/Track", + Gaudi::Property<std::string> m_trackOutputLoc + {this, "TracksOutputLocation", "Rec/Muon/Track", "LHCb tracks output location in TES"}; /// main steering reconstruction routine diff --git a/Muon/MuonTrackRec/src/component/MuonPadFromCoord.cpp b/Muon/MuonTrackRec/src/component/MuonPadFromCoord.cpp index 9ad18771522..9b6d93fb6ae 100644 --- a/Muon/MuonTrackRec/src/component/MuonPadFromCoord.cpp +++ b/Muon/MuonTrackRec/src/component/MuonPadFromCoord.cpp @@ -50,14 +50,34 @@ void MuonPadFromCoord::handle ( const Incident& incident ) } } - - - //============================================================================= StatusCode MuonPadFromCoord::initialize () { StatusCode sc = GaudiTool::initialize() ; if (!sc) return sc; + + // -- retrieve the pointer to the muon detector + if ( msgLevel(MSG::DEBUG) ) debug() << "Retrieve MuonDet" << endmsg; + m_muonDetector = getDet<DeMuonDetector>(DeMuonLocation::Default); + if(!m_muonDetector){ + error()<<"error retrieving Muon Detector geometry "<<endmsg; + return StatusCode::FAILURE; + } + + // geometry stuff ---> + auto nStations = m_muonDetector->stations(); + if ( nStations == 5 ) { + // Run 1 or 2 + m_stationOffset = 0; + } else if ( nStations == 4 ) { + // Upgrade + m_stationOffset = 1; + debug() << "Got " << nStations << " muon stations => setting offset to 1" << endmsg; + } else { + err() << "Wrong station number: " << nStations << endmsg; + return StatusCode::FAILURE; + } + incSvc()->addListener( this, IncidentType::EndEvent ); return StatusCode::SUCCESS ; } @@ -89,23 +109,23 @@ StatusCode MuonPadFromCoord::buildLogicalPads(const std::vector<MuonLogHit>& ) { m_pads.reserve( 2*coords->size() ); // loop over the coords for (const auto& coord : *coords) { - MuonTileID tile = MuonTileID(coord->key()); - unsigned int station = tile.station(); - unsigned int region = tile.region(); - double uncross = (station == 0 || ((station>2)&&(region==0))) ? false : coord->uncrossed(); + auto tile = MuonTileID(coord->key()); + auto station = tile.station() + m_stationOffset; // offset means 0 means M1 even in the upgrade + auto region = tile.region(); + auto uncross = (station == 0 || ((station>2)&&(region==0))) ? false : coord->uncrossed(); int time1= (int) coord->digitTDC1(); int time2= (int) coord->digitTDC2(); if(uncross) { - m_tiles.push_back(tile); + m_tiles.push_back( tile ); m_hits.emplace_back( tile ); - m_hits.back().setTime(time1); - m_pads.emplace_back(&m_hits.back()); // m_hits was reserved with enough capacity, so yes, this works + m_hits.back().setTime( time1 ); + m_pads.emplace_back( &m_hits.back() ); // m_hits was reserved with enough capacity, so yes, this works } else { const std::vector<MuonTileID > thetiles= coord->digitTile(); tile = MuonTileID(thetiles[0]); - m_tiles.push_back(tile); - m_hits.emplace_back(tile); + m_tiles.push_back( tile ); + m_hits.emplace_back( tile ); auto hit1 = &m_hits.back(); // m_hits was reserved with enough capacity, so yes, this works m_hits.back().setTime(time1); if (thetiles.size() == 1) { @@ -125,5 +145,3 @@ StatusCode MuonPadFromCoord::buildLogicalPads(const std::vector<MuonLogHit>& ) { m_padsReconstructed=true; return StatusCode::SUCCESS; } - - diff --git a/Muon/MuonTrackRec/src/component/MuonPadFromCoord.h b/Muon/MuonTrackRec/src/component/MuonPadFromCoord.h index 88f0923435a..a803cd3b96a 100644 --- a/Muon/MuonTrackRec/src/component/MuonPadFromCoord.h +++ b/Muon/MuonTrackRec/src/component/MuonPadFromCoord.h @@ -1,4 +1,4 @@ -#ifndef MUONPADFROMCOORD_H +#ifndef MUONPADFROMCOORD_H #define MUONPADFROMCOORD_H 1 #include <map> #include "GaudiAlg/GaudiTool.h" @@ -6,7 +6,7 @@ #include "MuonInterfaces/IMuonPadRec.h" // Interface /** @class MuonPadFromCoord MuonPadFromCoord.h - * converts MuonCoord to MuonLogPad objects + * converts MuonCoord to MuonLogPad objects * * @author G.Graziani * @date 2014-03-11 @@ -19,9 +19,9 @@ namespace LHCB { class MuonPadFromCoord final : public extends<GaudiTool, IMuonPadRec , IIncidentListener> { -public: +public: /// Standard constructor - MuonPadFromCoord( const std::string& type, + MuonPadFromCoord( const std::string& type, const std::string& name, const IInterface* parent); @@ -35,11 +35,13 @@ public: // from IIncidentListener void handle ( const Incident& incident ) override; private: - + std::vector<LHCb::MuonTileID> m_tiles; std::vector<MuonLogHit> m_hits; std::vector<MuonLogPad> m_pads; void clearPads(); + unsigned int m_stationOffset = 0; bool m_padsReconstructed = false; + DeMuonDetector* m_muonDetector = nullptr; }; #endif // MUONPADFROMCOORD_H diff --git a/Muon/MuonTrackRec/src/component/MuonTrackMomRec.cpp b/Muon/MuonTrackRec/src/component/MuonTrackMomRec.cpp index 83e78cf0e70..08cb6e45932 100644 --- a/Muon/MuonTrackRec/src/component/MuonTrackMomRec.cpp +++ b/Muon/MuonTrackRec/src/component/MuonTrackMomRec.cpp @@ -14,6 +14,7 @@ // from TrackEvent #include "Event/TrackParameters.h" +#include "Event/StateParameters.h" #include "Event/State.h" #include "Event/TrackTypes.h" #include "Event/Track.h" @@ -71,7 +72,7 @@ StatusCode MuonTrackMomRec::initBdl() { // compute integrated B field before M1 at x=y=0 // we will use this value for all tracks Gaudi::XYZPoint begin( 0., 0., 0. ); - Gaudi::XYZPoint end( 0., 0., m_muonDetector->getStationZ(0) ); + Gaudi::XYZPoint end( 0., 0., StateParameters::ZEndRich2 ); Gaudi::XYZVector bdl; m_bIntegrator -> calculateBdlAndCenter(begin, end, 0.0001, @@ -89,7 +90,7 @@ StatusCode MuonTrackMomRec::initBdl() { std::unique_ptr<LHCb::Track> MuonTrackMomRec::recMomentum(MuonTrack& track) const { StatusCode sc = StatusCode::SUCCESS; - double Zfirst = m_muonDetector->getStationZ(0); + const auto Zfirst = StateParameters::ZEndRich2; // create a state at the Z of M1 Gaudi::XYZPoint trackPos(track.bx() + track.sx() * Zfirst, track.by() + track.sy() * Zfirst, diff --git a/Tr/TrackTools/src/MuonSeeding.cpp b/Tr/TrackTools/src/MuonSeeding.cpp new file mode 100644 index 00000000000..f722b60be63 --- /dev/null +++ b/Tr/TrackTools/src/MuonSeeding.cpp @@ -0,0 +1,291 @@ +// Include files + +// local +#include "MuonSeeding.h" + +//----------------------------------------------------------------------------- +// Implementation file for class : MuonSeeding +// +// 2010-09-14 : Michel De Cian +//----------------------------------------------------------------------------- + +// Declaration of the Algorithm Factory +DECLARE_ALGORITHM_FACTORY( MuonSeeding ) + +//============================================================================= +// Initialization +//============================================================================= +StatusCode MuonSeeding::initialize() { + const StatusCode sc = GaudiAlgorithm::initialize(); // must be executed first + if ( sc.isFailure() ) { + return sc; // error printed already by GaudiAlgorithm + } + + if ( msgLevel(MSG::DEBUG) ) { + debug() << "==> Initialize" << endmsg; + } + + // -- Tracking tools + m_trackTool = tool<IMuonTrackRec>( m_trackToolName, this ); + m_extrapolator = tool<ITrackExtrapolator>( m_extrapolatorName, this ); + m_trackFitter = tool<ITrackFitter>("TrackMasterFitter",this); + m_momentumTool = tool<IMuonTrackMomRec>("MuonTrackMomRec", this); + + // -- MC tools + if(m_MC){ + m_muonPad2MC = tool<IMuonPad2MCTool>("MuonPad2MCTool", this); + m_lhcbid2mcparticles = tool<ILHCbIDsToMCParticles>("LHCbIDsToMCParticles", this); + } + + return sc; +} + +//============================================================================= +// Main execution +//============================================================================= +StatusCode MuonSeeding::execute() +{ + //Chrono chrono( chronoSvc(),name()+"::execute()" ); + + if ( msgLevel(MSG::DEBUG) ) debug() << "==> Execute" << endmsg; + + // -- This is where we assume the track came from + std::vector<double> PVPos = { 0.0, 0.0, 0.0 }; + + const std::vector<MuonTrack>& muonTracks = m_trackTool->tracks(); + auto tracks = std::make_unique<LHCb::Tracks>(); + + // -- Loop over all Muon Tracks + for (auto it = muonTracks.begin() ; it != muonTracks.end() ; it++){ + + const LHCb::MCParticle* bestMCPart = nullptr; + + // -- Associate Muon Hits to an MCParticle, get momentum components + double truePx = 0; + double truePy = 0; + double truePz = 0; + double trueP = 0; + + if(m_MC){ + bestMCPart = assocMCParticle(it->getHits()); + if(bestMCPart){ + truePx = bestMCPart->momentum().X(); + truePy = bestMCPart->momentum().Y(); + truePz = bestMCPart->momentum().Z(); + trueP = bestMCPart->p(); + } + } + // --------------------------------------------------------------- + + + // -- Get the momentum of the muon track, make a new track + auto track = m_momentumTool->recMomentum( const_cast<MuonTrack&>(*it) ); //@FIXME + if(!track) continue; + // --------------------------------------------------------------- + + // -- Change Q/p until it points to the PV (stolen from Wouter) + LHCb::State muonState; + auto sc = iterateToPV(track.get(), muonState, PVPos, it->qOverP() ); // -- This is the function that iterates + if(!sc){ + Warning( "==> Could not iterate to PV!", StatusCode::SUCCESS, 0 ).ignore(); + continue; + } + + // -- Set Pattern Reco status and track type, finally fit the track + track->setPatRecStatus( LHCb::Track::PatRecStatus::PatRecIDs ); + track->setType( LHCb::Track::Types::Muon ); + track->clearStates( ); // remove the state recMomentum created + track->addToStates( muonState ); + // --------------------------------------------------------------- + + // -- Add some MC info to the track + if(m_MC){ + track->addInfo(1300, truePx); + track->addInfo(1301, truePy); + track->addInfo(1302, truePz); + track->addInfo(1303, trueP); + } + + // -- Finally, insert the track into the container! + tracks->insert( track.release() ); + } + + // --------------------------------------------------------------- + + if(msgLevel(MSG::DEBUG)) { + debug() << "Filling " << tracks->size() << " tracks in " << m_outputLoc << endmsg; + } + + bool filterPassed = false; + + if ( !tracks->empty() ) { + put( tracks.release(), m_outputLoc ); + filterPassed = true; + } + + // -- Only go on if MuonTT tracks have been reconstructed + setFilterPassed( filterPassed ); + + return StatusCode::SUCCESS; +} + +//============================================================================= +// Fill Muon Stub Information +//============================================================================= +void MuonSeeding::fillMuonStubInfo(LHCb::Track& track, const MuonTrack& muTrack) const{ + + // -- Retrieve fit errors + Gaudi::XYZVector slope; + Gaudi::SymMatrix3x3 errmat; + track.slopes( slope , errmat ); + + + // -- Get parameters of the MuonStub + double bx_fit = track.position().X() - slope.X()*track.position().Z(); + double by_fit = track.position().Y() - slope.Y()*track.position().Z(); + + double sx_fit = slope.X(); + double sy_fit = slope.Y(); + + double chi2_fit = track.chi2PerDoF(); + + double npad_fit = double( muTrack.getHits().size() ); + // --------------------------------------------------------------- + + double sx = muTrack.sx() ; // -- track slope in XZ + double bx = muTrack.bx() ; // -- track intercept in XZ + double sy = muTrack.sy() ; // -- track slope in YZ + double by = muTrack.by() ; // -- track intercept in YZ + double esx = muTrack.errsx(); // -- error on sx + double ebx = muTrack.errbx(); // -- error on bx + double esy = muTrack.errsy(); // -- error on sy + double eby = muTrack.errby(); // -- error on by + + // --------------------------------------------------------------- + track.addInfo(1000, muTrack.chi2x() ); + track.addInfo(1001, muTrack.chi2y() ); + track.addInfo(1002, sx); + track.addInfo(1003, bx); + track.addInfo(1004, sy); + track.addInfo(1005, by); + track.addInfo(1006, esx); + track.addInfo(1007, ebx); + track.addInfo(1008, esy); + track.addInfo(1009, eby); + track.addInfo(1010, bx_fit); + track.addInfo(1011, by_fit); + track.addInfo(1012, sx_fit); + track.addInfo(1013, sy_fit); + track.addInfo(1014, chi2_fit); + track.addInfo(1015, npad_fit); + track.addInfo(1016, muTrack.p()); + track.addInfo(1017, muTrack.pt()); + track.addInfo(1018, muTrack.qOverP()); + track.addInfo(1019, muTrack.covbsy()); + + track.addInfo(1020, muTrack.getSpan() ); + track.addInfo(1021, muTrack.qOverP() ); + + track.addInfo(1200, muTrack.momentum().X()); + track.addInfo(1201, muTrack.momentum().Y()); + track.addInfo(1202, muTrack.momentum().Z()); + // --------------------------------------------------------------- + +} +//============================================================================= +// Change the q/p till the track points to the PV (stolen from Wouter) +//============================================================================= +StatusCode MuonSeeding::iterateToPV(LHCb::Track* track, LHCb::State& finalState, const std::vector + double> &PVPos, double qOverP) +{ + LHCb::State muonState = track->closestState( 15000 ); + + double dXdQoP = 1e7 ; + double deltaX = 100; + muonState.setQOverP( qOverP ) ; + + // -- Now call the extrapolator and iterate until we have the desired accuracy. + double tolerance = 0.5 ; // [mm] + Gaudi::TrackMatrix jacobian ; + + LHCb::State dummyState = muonState; + + for(int i=0; i<10 && std::abs(deltaX) > tolerance ; ++i) { + dummyState = muonState ; + StatusCode sc = m_extrapolator->propagate( dummyState, PVPos[2], &jacobian) ; + if(!sc) return StatusCode::FAILURE; + dXdQoP = jacobian(0,4) ; + deltaX = -(dummyState.x() - PVPos[0]) ; + double deltaQoP = deltaX / dXdQoP ; + muonState.setQOverP( muonState.qOverP() + deltaQoP ) ; + } + + finalState = muonState; + return StatusCode::SUCCESS; + +} +//============================================================================= +// Associate MCParticle to Pad Hits +//============================================================================= +const LHCb::MCParticle* MuonSeeding::assocMCParticle(std::vector<const MuonHit*> muonHits){ + + std::map<const LHCb::MCParticle*, int> mcparts; + + // -- Get the MCParticles which are associated to the hits of the Particle in the Muon System + // -- Put the mcparticles into a map, count how many times they appear + for( const auto& muonHit : muonHits ) { + const LHCb::MCParticle *p = m_muonPad2MC->Pad2MC( muonHit->centerTile() ); + if (!p) continue; + auto iter = mcparts.find(p); + if( iter != mcparts.end()){ + ++(iter->second); + } else { + mcparts.emplace( p, 1 ); + } + } + // --------------------------------------------------------------- + + // -- Get particle which contributes the most (i.e. appears the most) + std::pair<const LHCb::MCParticle*, int> best{nullptr, 0}; + for(const auto& contrib : mcparts) { + if( contrib.second > best.second ){ + best.first = contrib.first; + best.second = contrib.second; + } + } + // --------------------------------------------------------------- + + // -- Call a Particle associated if more than 70% of all MCParticles for all the hits on a track are the same + bool isTrue = false; + if( double(best.second) / muonHits.size() > 0.7){ + + // -- Particle has to be a Muon + if( best.first->particleID().abspid() == 13 ){ + isTrue = true; + } + }else{ + if(msgLevel(MSG::DEBUG)) debug() << "could not assign particle" << endmsg; + } + // --------------------------------------------------------------- + + return isTrue ? best.first : nullptr; + +} +//========================================================================= +// Check if a given LHCbID is on "track" of a MCParticle (not called anywhere at the moment, +// for occasional testing needed) +//========================================================================= +bool MuonSeeding::isIDOnTrack(LHCb::LHCbID id, LHCb::MCParticle* mcp){ + + + std::map<LHCb::MCParticle*, unsigned int> linkMap; + // -- This class should link my id to the MCParticle! + m_lhcbid2mcparticles->link(id, linkMap); + + // -- For simplicity, we only look at unambigously assigned mcparticles! + if(linkMap.size() > 1 ) return false; + + return std::any_of( linkMap.begin(), linkMap.end(), + [&](const std::pair<const LHCb::MCParticle*,unsigned int>& i) + { return mcp==i.first; }); +} diff --git a/Tr/TrackTools/src/MuonSeeding.h b/Tr/TrackTools/src/MuonSeeding.h new file mode 100644 index 00000000000..507045e5dfb --- /dev/null +++ b/Tr/TrackTools/src/MuonSeeding.h @@ -0,0 +1,76 @@ +#ifndef MUONSEEDING_H +#define MUONSEEDING_H 1 + +// Include files +// from Gaudi +#include "GaudiAlg/GaudiAlgorithm.h" +#include "MuonInterfaces/MuonTrack.h" +#include "MuonInterfaces/MuonHit.h" +#include "MuonInterfaces/IMuonTrackMomRec.h" +#include "MuonInterfaces/IMuonTrackRec.h" +#include "Event/Track.h" +#include "Event/Particle.h" + +#include "MCInterfaces/IMuonPad2MCTool.h" +#include "Event/MCParticle.h" +#include "Event/State.h" +#include "Event/RecVertex.h" + +#include "TrackInterfaces/ITrackFitter.h" +#include "TrackInterfaces/ITrackExtrapolator.h" + +#include "Kernel/LHCbID.h" +#include "MCInterfaces/ILHCbIDsToMCParticles.h" + +#include "GaudiKernel/Chrono.h" + + +/** @class MuonSeeding MuonSeeding.h + * + * \brief Make a MuonSeeding: Get muon standalone tracks + * + * Parameters: + * - ToolName: Name for the tool that makes muon standalone track. + * - Extrapolator: Name for the track extrapolator. + * - MC: To enable MC association. + * - FillMuonStubInfo: Fill parameters of muon stub in info fields of track; + * - Output: The location the tracks should be written to. + * + * @author Michel De Cian + * @date 2010-09-20 + */ + +class MuonSeeding final : public GaudiAlgorithm { +public: + /// Standard constructor + using GaudiAlgorithm::GaudiAlgorithm; + + StatusCode initialize() override; ///< Algorithm initialization + StatusCode execute() override; ///< Algorithm execution + +private: + + // -- Methods + StatusCode fillPVs(std::vector<double>& PVPos); + void fillMuonStubInfo(LHCb::Track& track, const MuonTrack& muTrack) const; + StatusCode iterateToPV(LHCb::Track* track, LHCb::State& finalState, std::vector<double> PVPos, double qOverP); + const LHCb::MCParticle* assocMCParticle(const std::vector<const MuonHit*> muonHits); + bool isIDOnTrack(LHCb::LHCbID id, LHCb::MCParticle* mcp); + + // -- Properties + Gaudi::Property<std::string> m_outputLoc { this, "Output", "Rec/"+name()+"/Tracks" }; + Gaudi::Property<std::string> m_pvLoc { this, "PVLocation", LHCb::RecVertexLocation::Primary }; + Gaudi::Property<std::string> m_extrapolatorName { this, "Extrapolator", "TrackMasterExtrapolator" }; + Gaudi::Property<bool> m_MC { this, "MC", false }; + Gaudi::Property<bool> m_fillMuonStubInfo { this, "FillMuonStubInfo", false }; + Gaudi::Property<std::string> m_trackToolName { this, "ToolName", "MuonNNetRec" }; + + // -- Tools + ILHCbIDsToMCParticles* m_lhcbid2mcparticles = nullptr; + IMuonPad2MCTool* m_muonPad2MC = nullptr; + IMuonTrackRec* m_trackTool = nullptr; + IMuonTrackMomRec* m_momentumTool = nullptr; + ITrackExtrapolator* m_extrapolator = nullptr; + ITrackFitter* m_trackFitter = nullptr; +}; +#endif // GETMUONTRACK_H -- GitLab From 1677f2cab133e5a5d1965552eff4e8121353867e Mon Sep 17 00:00:00 2001 From: Olli Lupton <oliver.lupton@cern.ch> Date: Tue, 14 Nov 2017 23:49:46 +0100 Subject: [PATCH 2/5] Rejig MC matching/linking for the muon stubs. --- Pr/PrMCTools/src/PrMuonStubAssociator.cpp | 108 ++++++++++++++++++++ Pr/PrMCTools/src/PrMuonStubAssociator.h | 41 ++++++++ Tr/TrackTools/src/MuonSeeding.cpp | 116 +--------------------- Tr/TrackTools/src/MuonSeeding.h | 12 +-- Tr/TrackTools/src/TrackInitFit.h | 2 +- 5 files changed, 156 insertions(+), 123 deletions(-) create mode 100644 Pr/PrMCTools/src/PrMuonStubAssociator.cpp create mode 100644 Pr/PrMCTools/src/PrMuonStubAssociator.h diff --git a/Pr/PrMCTools/src/PrMuonStubAssociator.cpp b/Pr/PrMCTools/src/PrMuonStubAssociator.cpp new file mode 100644 index 00000000000..59ba6f8dbe4 --- /dev/null +++ b/Pr/PrMCTools/src/PrMuonStubAssociator.cpp @@ -0,0 +1,108 @@ +// Include files + + +// from Gaudi +#include "GaudiKernel/IDataManagerSvc.h" +#include "GaudiKernel/SmartIF.h" + +// from Event +#include "Event/Track.h" +#include "Linker/LinkerWithKey.h" + +#include "Event/MCParticle.h" +#include "Event/MCVertex.h" + +// local +#include "PrMuonStubAssociator.h" + +DECLARE_ALGORITHM_FACTORY( PrMuonStubAssociator ) + +//============================================================================= +// Standard constructor, initializes variables +//============================================================================= + PrMuonStubAssociator::PrMuonStubAssociator( const std::string& name, ISvcLocator* pSvcLocator ) + : GaudiAlgorithm( name, pSvcLocator ) +{ +} + +//============================================================================= +// Destructor +//============================================================================= +PrMuonStubAssociator::~PrMuonStubAssociator() { } + +//============================================================================= +// Initialisation. Check parameters +//============================================================================= +StatusCode PrMuonStubAssociator::initialize() +{ + // Mandatory initialization of GaudiAlgorithm + StatusCode sc = GaudiAlgorithm::initialize(); + if( sc.isFailure() ) { + return sc; + } + + info() << "Processing container: " << m_container.value() << endmsg; + + m_muonPad2MC = tool<IMuonPad2MCTool>("MuonPad2MCTool", this); + if ( !m_muonPad2MC ) { + error() << "Failed to create muonPad2MC tool!" << endmsg; + return StatusCode::FAILURE; + } + + return StatusCode::SUCCESS; +} + +//============================================================================= +// Main execution +//============================================================================= +StatusCode PrMuonStubAssociator::execute() +{ + // Retrieve the MCParticles + auto mcParts = getIfExists<LHCb::MCParticles> ( LHCb::MCParticleLocation::Default ); + if ( !mcParts ) { + if( msgLevel(MSG::ERROR) ) { + error() << "MCParticles at " << LHCb::MCParticleLocation::Default << "' do not exist" <<endmsg; + } + return StatusCode::SUCCESS; + } + + // Create the Linker table from Track to MCParticle + // Sorted by decreasing weight, so first retrieved has highest weight + // This has to be done, even if there are no tracks in the event, to satisfy the DST writer + LinkerWithKey<LHCb::MCParticle,LHCb::Track> myLinker( evtSvc(), msgSvc(), m_container ); + myLinker.reset(); + + // Retrieve the Tracks + auto tracks = getIfExists<LHCb::Track::Range>( m_container ); + if ( tracks.empty() ) { + if( msgLevel(MSG::DEBUG) ) { + debug() << "No tracks in container '" << m_container.value() << "'. Skipping." <<endmsg; + } + return StatusCode::SUCCESS; + } else if( msgLevel(MSG::DEBUG) ) { + debug() << "Loaded " << tracks.size() << " tracks from '" << m_container.value() << "'" << endmsg; + } + + for ( const auto& track : tracks ) { + auto nMuon = 0; + std::map<LHCb::MCParticle*,unsigned int> matches; + for ( const auto& id : track->lhcbIDs() ) { + if ( id.isMuon() ) { + nMuon++; + auto mc_part = m_muonPad2MC->Pad2MC( id.muonID() ); + if ( mc_part && mc_part->parent() == mcParts ) { + matches[ mc_part ]++; + } + } + } + + for ( const auto& match : matches ) { + auto fraction = 1.0 * match.second / nMuon; + if ( fraction > m_matchThreshold ) { + myLinker.link( track, match.first, fraction ); + } + } + } + + return StatusCode::SUCCESS; +} diff --git a/Pr/PrMCTools/src/PrMuonStubAssociator.h b/Pr/PrMCTools/src/PrMuonStubAssociator.h new file mode 100644 index 00000000000..8ac429d6a73 --- /dev/null +++ b/Pr/PrMCTools/src/PrMuonStubAssociator.h @@ -0,0 +1,41 @@ +#ifndef PRMUONSTUBASSOCIATOR_H +#define PRMUONSTUBASSOCIATOR_H 1 + +// Include files +#include "GaudiAlg/GaudiAlgorithm.h" +#include "MCInterfaces/IMuonPad2MCTool.h" + +/** @class TrackAssociator TrackAssociator.h + * + * This algorithm computes the link between a Track and a MCParticle. + * The requirement is a match of both the Velo/VP and the T part of the + * Track. If there are not enough coordinates, the match is assumed so that + * a Velo only or a T only are matched properly. + * The required fraction of hits is a jobOption 'FractionOK', default 0.70. + * + * Rewritten for the upgrade, handles all containers in one instance + * + * @author Olivier Callot + * @date 2012-04-04 + */ + +class PrMuonStubAssociator : public GaudiAlgorithm { +public: + + // Standard constructor + PrMuonStubAssociator( const std::string& name, ISvcLocator* pSvcLocator ); + + // Destructor + virtual ~PrMuonStubAssociator(); + + StatusCode initialize() override; ///< Algorithm initialization + StatusCode execute() override; ///< Algorithm execution + +private: + Gaudi::Property<std::string> m_container { this, "Container", "Rec/Track/StandaloneMuonFitted" }; + Gaudi::Property<float> m_matchThreshold { this, "MatchThreshold", 0.7 }; + + IMuonPad2MCTool* m_muonPad2MC = nullptr; +}; + +#endif // PRMUONSTUBASSOCIATOR_H diff --git a/Tr/TrackTools/src/MuonSeeding.cpp b/Tr/TrackTools/src/MuonSeeding.cpp index f722b60be63..5c2d53185a1 100644 --- a/Tr/TrackTools/src/MuonSeeding.cpp +++ b/Tr/TrackTools/src/MuonSeeding.cpp @@ -31,12 +31,6 @@ StatusCode MuonSeeding::initialize() { m_trackFitter = tool<ITrackFitter>("TrackMasterFitter",this); m_momentumTool = tool<IMuonTrackMomRec>("MuonTrackMomRec", this); - // -- MC tools - if(m_MC){ - m_muonPad2MC = tool<IMuonPad2MCTool>("MuonPad2MCTool", this); - m_lhcbid2mcparticles = tool<ILHCbIDsToMCParticles>("LHCbIDsToMCParticles", this); - } - return sc; } @@ -57,27 +51,6 @@ StatusCode MuonSeeding::execute() // -- Loop over all Muon Tracks for (auto it = muonTracks.begin() ; it != muonTracks.end() ; it++){ - - const LHCb::MCParticle* bestMCPart = nullptr; - - // -- Associate Muon Hits to an MCParticle, get momentum components - double truePx = 0; - double truePy = 0; - double truePz = 0; - double trueP = 0; - - if(m_MC){ - bestMCPart = assocMCParticle(it->getHits()); - if(bestMCPart){ - truePx = bestMCPart->momentum().X(); - truePy = bestMCPart->momentum().Y(); - truePz = bestMCPart->momentum().Z(); - trueP = bestMCPart->p(); - } - } - // --------------------------------------------------------------- - - // -- Get the momentum of the muon track, make a new track auto track = m_momentumTool->recMomentum( const_cast<MuonTrack&>(*it) ); //@FIXME if(!track) continue; @@ -98,14 +71,6 @@ StatusCode MuonSeeding::execute() track->addToStates( muonState ); // --------------------------------------------------------------- - // -- Add some MC info to the track - if(m_MC){ - track->addInfo(1300, truePx); - track->addInfo(1301, truePy); - track->addInfo(1302, truePz); - track->addInfo(1303, trueP); - } - // -- Finally, insert the track into the container! tracks->insert( track.release() ); } @@ -116,15 +81,11 @@ StatusCode MuonSeeding::execute() debug() << "Filling " << tracks->size() << " tracks in " << m_outputLoc << endmsg; } - bool filterPassed = false; + // Control flow can stop if we didn't make anything + setFilterPassed( !tracks->empty() ); - if ( !tracks->empty() ) { - put( tracks.release(), m_outputLoc ); - filterPassed = true; - } - - // -- Only go on if MuonTT tracks have been reconstructed - setFilterPassed( filterPassed ); + // Put our output on the TES + put( tracks.release(), m_outputLoc ); return StatusCode::SUCCESS; } @@ -195,8 +156,7 @@ void MuonSeeding::fillMuonStubInfo(LHCb::Track& track, const MuonTrack& muTrack) //============================================================================= // Change the q/p till the track points to the PV (stolen from Wouter) //============================================================================= -StatusCode MuonSeeding::iterateToPV(LHCb::Track* track, LHCb::State& finalState, const std::vector - double> &PVPos, double qOverP) +StatusCode MuonSeeding::iterateToPV(LHCb::Track* track, LHCb::State& finalState, const std::vector<double> &PVPos, double qOverP) { LHCb::State muonState = track->closestState( 15000 ); @@ -222,70 +182,4 @@ StatusCode MuonSeeding::iterateToPV(LHCb::Track* track, LHCb::State& finalState, finalState = muonState; return StatusCode::SUCCESS; - -} -//============================================================================= -// Associate MCParticle to Pad Hits -//============================================================================= -const LHCb::MCParticle* MuonSeeding::assocMCParticle(std::vector<const MuonHit*> muonHits){ - - std::map<const LHCb::MCParticle*, int> mcparts; - - // -- Get the MCParticles which are associated to the hits of the Particle in the Muon System - // -- Put the mcparticles into a map, count how many times they appear - for( const auto& muonHit : muonHits ) { - const LHCb::MCParticle *p = m_muonPad2MC->Pad2MC( muonHit->centerTile() ); - if (!p) continue; - auto iter = mcparts.find(p); - if( iter != mcparts.end()){ - ++(iter->second); - } else { - mcparts.emplace( p, 1 ); - } - } - // --------------------------------------------------------------- - - // -- Get particle which contributes the most (i.e. appears the most) - std::pair<const LHCb::MCParticle*, int> best{nullptr, 0}; - for(const auto& contrib : mcparts) { - if( contrib.second > best.second ){ - best.first = contrib.first; - best.second = contrib.second; - } - } - // --------------------------------------------------------------- - - // -- Call a Particle associated if more than 70% of all MCParticles for all the hits on a track are the same - bool isTrue = false; - if( double(best.second) / muonHits.size() > 0.7){ - - // -- Particle has to be a Muon - if( best.first->particleID().abspid() == 13 ){ - isTrue = true; - } - }else{ - if(msgLevel(MSG::DEBUG)) debug() << "could not assign particle" << endmsg; - } - // --------------------------------------------------------------- - - return isTrue ? best.first : nullptr; - -} -//========================================================================= -// Check if a given LHCbID is on "track" of a MCParticle (not called anywhere at the moment, -// for occasional testing needed) -//========================================================================= -bool MuonSeeding::isIDOnTrack(LHCb::LHCbID id, LHCb::MCParticle* mcp){ - - - std::map<LHCb::MCParticle*, unsigned int> linkMap; - // -- This class should link my id to the MCParticle! - m_lhcbid2mcparticles->link(id, linkMap); - - // -- For simplicity, we only look at unambigously assigned mcparticles! - if(linkMap.size() > 1 ) return false; - - return std::any_of( linkMap.begin(), linkMap.end(), - [&](const std::pair<const LHCb::MCParticle*,unsigned int>& i) - { return mcp==i.first; }); } diff --git a/Tr/TrackTools/src/MuonSeeding.h b/Tr/TrackTools/src/MuonSeeding.h index 507045e5dfb..394c545444e 100644 --- a/Tr/TrackTools/src/MuonSeeding.h +++ b/Tr/TrackTools/src/MuonSeeding.h @@ -11,8 +11,6 @@ #include "Event/Track.h" #include "Event/Particle.h" -#include "MCInterfaces/IMuonPad2MCTool.h" -#include "Event/MCParticle.h" #include "Event/State.h" #include "Event/RecVertex.h" @@ -20,8 +18,6 @@ #include "TrackInterfaces/ITrackExtrapolator.h" #include "Kernel/LHCbID.h" -#include "MCInterfaces/ILHCbIDsToMCParticles.h" - #include "GaudiKernel/Chrono.h" @@ -32,7 +28,6 @@ * Parameters: * - ToolName: Name for the tool that makes muon standalone track. * - Extrapolator: Name for the track extrapolator. - * - MC: To enable MC association. * - FillMuonStubInfo: Fill parameters of muon stub in info fields of track; * - Output: The location the tracks should be written to. * @@ -53,21 +48,16 @@ private: // -- Methods StatusCode fillPVs(std::vector<double>& PVPos); void fillMuonStubInfo(LHCb::Track& track, const MuonTrack& muTrack) const; - StatusCode iterateToPV(LHCb::Track* track, LHCb::State& finalState, std::vector<double> PVPos, double qOverP); - const LHCb::MCParticle* assocMCParticle(const std::vector<const MuonHit*> muonHits); - bool isIDOnTrack(LHCb::LHCbID id, LHCb::MCParticle* mcp); + StatusCode iterateToPV(LHCb::Track* track, LHCb::State& finalState, const std::vector<double>& PVPos, double qOverP); // -- Properties Gaudi::Property<std::string> m_outputLoc { this, "Output", "Rec/"+name()+"/Tracks" }; Gaudi::Property<std::string> m_pvLoc { this, "PVLocation", LHCb::RecVertexLocation::Primary }; Gaudi::Property<std::string> m_extrapolatorName { this, "Extrapolator", "TrackMasterExtrapolator" }; - Gaudi::Property<bool> m_MC { this, "MC", false }; Gaudi::Property<bool> m_fillMuonStubInfo { this, "FillMuonStubInfo", false }; Gaudi::Property<std::string> m_trackToolName { this, "ToolName", "MuonNNetRec" }; // -- Tools - ILHCbIDsToMCParticles* m_lhcbid2mcparticles = nullptr; - IMuonPad2MCTool* m_muonPad2MC = nullptr; IMuonTrackRec* m_trackTool = nullptr; IMuonTrackMomRec* m_momentumTool = nullptr; ITrackExtrapolator* m_extrapolator = nullptr; diff --git a/Tr/TrackTools/src/TrackInitFit.h b/Tr/TrackTools/src/TrackInitFit.h index dd66d9b3c00..e47cc259e65 100644 --- a/Tr/TrackTools/src/TrackInitFit.h +++ b/Tr/TrackTools/src/TrackInitFit.h @@ -33,7 +33,7 @@ public: // and fit return m_fitTrack->operator()(track, pid); } - + StatusCode operator() ( std::vector<std::reference_wrapper<LHCb::Track>>& tracks, const LHCb::Tr::PID& pid -- GitLab From d3db47d04cdd0346885240c4378fdaf50ac475b6 Mon Sep 17 00:00:00 2001 From: Olli Lupton <oliver.lupton@cern.ch> Date: Wed, 15 Nov 2017 13:47:58 +0100 Subject: [PATCH 3/5] Add an option to fit the muon stubs before massaging them to point at the origin. --- Tr/TrackTools/src/MuonSeeding.cpp | 175 +++++++++++++++--------------- Tr/TrackTools/src/MuonSeeding.h | 20 ++-- 2 files changed, 96 insertions(+), 99 deletions(-) diff --git a/Tr/TrackTools/src/MuonSeeding.cpp b/Tr/TrackTools/src/MuonSeeding.cpp index 5c2d53185a1..4024a48241a 100644 --- a/Tr/TrackTools/src/MuonSeeding.cpp +++ b/Tr/TrackTools/src/MuonSeeding.cpp @@ -2,6 +2,7 @@ // local #include "MuonSeeding.h" +#include "Event/VPCluster.h" //----------------------------------------------------------------------------- // Implementation file for class : MuonSeeding @@ -12,6 +13,15 @@ // Declaration of the Algorithm Factory DECLARE_ALGORITHM_FACTORY( MuonSeeding ) +MuonSeeding::MuonSeeding( const std::string& name, ISvcLocator* pSvcLocator ) + : GaudiAlgorithm( name, pSvcLocator ) +{ + declareProperty( "MuonRecTool", m_trackTool ); + declareProperty( "Fitter", m_trackFitter ); + declareProperty( "MuomMomTool", m_momentumTool ); + declareProperty( "Extrapolator", m_extrapolator ); +} + //============================================================================= // Initialization //============================================================================= @@ -26,10 +36,10 @@ StatusCode MuonSeeding::initialize() { } // -- Tracking tools - m_trackTool = tool<IMuonTrackRec>( m_trackToolName, this ); - m_extrapolator = tool<ITrackExtrapolator>( m_extrapolatorName, this ); - m_trackFitter = tool<ITrackFitter>("TrackMasterFitter",this); - m_momentumTool = tool<IMuonTrackMomRec>("MuonTrackMomRec", this); + m_trackTool.retrieve(); + m_extrapolator.retrieve(); + m_trackFitter.retrieve(); + m_momentumTool.retrieve(); return sc; } @@ -44,22 +54,58 @@ StatusCode MuonSeeding::execute() if ( msgLevel(MSG::DEBUG) ) debug() << "==> Execute" << endmsg; // -- This is where we assume the track came from - std::vector<double> PVPos = { 0.0, 0.0, 0.0 }; + std::vector<double> PVPos = { 0.0, 0.0, 0.0 }; // TODO take the VELO position? + + /*const LHCb::VPCluster* vp_cluster = nullptr; + if ( true ) { // VP cluster 'hack' + auto vp_clusters = getIfExists<LHCb::VPClusters>( LHCb::VPClusterLocation::Default ); + if ( !vp_clusters || vp_clusters->empty() ) { + if ( msgLevel(MSG::DEBUG) ) { + debug() << "No VP clusters, skipping..." << endmsg; + } + + return StatusCode::SUCCESS; + } + + // Find the VP cluster nearest the origin + std::vector<std::pair<float,const LHCb::VPCluster*> > scores; + for ( const auto& cluster : *vp_clusters ) { + scores.emplace_back( + std::sqrt( std::pow( cluster->x(), 2) + std::pow( cluster->y(), 2) + std::pow( cluster->z(), 2) ), + cluster ); + } - const std::vector<MuonTrack>& muonTracks = m_trackTool->tracks(); + auto min_iter = std::min_element( scores.begin(), scores.end() ); + vp_cluster = min_iter->second; + }*/ + + const auto& muonTracks = m_trackTool->tracks(); auto tracks = std::make_unique<LHCb::Tracks>(); // -- Loop over all Muon Tracks - for (auto it = muonTracks.begin() ; it != muonTracks.end() ; it++){ + for ( const auto& muonTrack : muonTracks ) { // -- Get the momentum of the muon track, make a new track - auto track = m_momentumTool->recMomentum( const_cast<MuonTrack&>(*it) ); //@FIXME - if(!track) continue; - // --------------------------------------------------------------- + auto track = m_momentumTool->recMomentum( const_cast<MuonTrack&>( muonTrack ) ); //@FIXME + if ( !track ) { + continue; + } + + if ( m_fitTracks ) { + // Try and improve the covariance information + auto sc = m_trackFitter->fit( *track ); + if ( sc.isFailure() ) { + if ( msgLevel( MSG::WARNING ) ) { + warning() << "Track fit failed" << endmsg; + } + + continue; + } + } - // -- Change Q/p until it points to the PV (stolen from Wouter) - LHCb::State muonState; - auto sc = iterateToPV(track.get(), muonState, PVPos, it->qOverP() ); // -- This is the function that iterates - if(!sc){ + // -- Change q/p until it points to the origin (adapted from Wouter) + LHCb::State veloState, muonState; + auto sc = iterateToPV(track.get(), muonState, veloState, PVPos, muonTrack.qOverP() ); + if ( sc.isFailure() ) { Warning( "==> Could not iterate to PV!", StatusCode::SUCCESS, 0 ).ignore(); continue; } @@ -67,7 +113,12 @@ StatusCode MuonSeeding::execute() // -- Set Pattern Reco status and track type, finally fit the track track->setPatRecStatus( LHCb::Track::PatRecStatus::PatRecIDs ); track->setType( LHCb::Track::Types::Muon ); - track->clearStates( ); // remove the state recMomentum created + + if ( !m_fitTracks ) { + track->clearStates( ); // remove the state recMomentum created + } + + track->addToStates( veloState ); track->addToStates( muonState ); // --------------------------------------------------------------- @@ -90,96 +141,44 @@ StatusCode MuonSeeding::execute() return StatusCode::SUCCESS; } -//============================================================================= -// Fill Muon Stub Information -//============================================================================= -void MuonSeeding::fillMuonStubInfo(LHCb::Track& track, const MuonTrack& muTrack) const{ - - // -- Retrieve fit errors - Gaudi::XYZVector slope; - Gaudi::SymMatrix3x3 errmat; - track.slopes( slope , errmat ); - - - // -- Get parameters of the MuonStub - double bx_fit = track.position().X() - slope.X()*track.position().Z(); - double by_fit = track.position().Y() - slope.Y()*track.position().Z(); - - double sx_fit = slope.X(); - double sy_fit = slope.Y(); - - double chi2_fit = track.chi2PerDoF(); - - double npad_fit = double( muTrack.getHits().size() ); - // --------------------------------------------------------------- - - double sx = muTrack.sx() ; // -- track slope in XZ - double bx = muTrack.bx() ; // -- track intercept in XZ - double sy = muTrack.sy() ; // -- track slope in YZ - double by = muTrack.by() ; // -- track intercept in YZ - double esx = muTrack.errsx(); // -- error on sx - double ebx = muTrack.errbx(); // -- error on bx - double esy = muTrack.errsy(); // -- error on sy - double eby = muTrack.errby(); // -- error on by - - // --------------------------------------------------------------- - track.addInfo(1000, muTrack.chi2x() ); - track.addInfo(1001, muTrack.chi2y() ); - track.addInfo(1002, sx); - track.addInfo(1003, bx); - track.addInfo(1004, sy); - track.addInfo(1005, by); - track.addInfo(1006, esx); - track.addInfo(1007, ebx); - track.addInfo(1008, esy); - track.addInfo(1009, eby); - track.addInfo(1010, bx_fit); - track.addInfo(1011, by_fit); - track.addInfo(1012, sx_fit); - track.addInfo(1013, sy_fit); - track.addInfo(1014, chi2_fit); - track.addInfo(1015, npad_fit); - track.addInfo(1016, muTrack.p()); - track.addInfo(1017, muTrack.pt()); - track.addInfo(1018, muTrack.qOverP()); - track.addInfo(1019, muTrack.covbsy()); - - track.addInfo(1020, muTrack.getSpan() ); - track.addInfo(1021, muTrack.qOverP() ); - - track.addInfo(1200, muTrack.momentum().X()); - track.addInfo(1201, muTrack.momentum().Y()); - track.addInfo(1202, muTrack.momentum().Z()); - // --------------------------------------------------------------- - -} //============================================================================= // Change the q/p till the track points to the PV (stolen from Wouter) //============================================================================= -StatusCode MuonSeeding::iterateToPV(LHCb::Track* track, LHCb::State& finalState, const std::vector<double> &PVPos, double qOverP) +StatusCode MuonSeeding::iterateToPV(LHCb::Track* track, LHCb::State& muonState, LHCb::State& veloState, + const std::vector<double> &PVPos, double qOverP) { - LHCb::State muonState = track->closestState( 15000 ); + muonState = track->closestState( 15000 ); double dXdQoP = 1e7 ; double deltaX = 100; muonState.setQOverP( qOverP ) ; + // Set the y slope based on the target position at ~the origin + muonState.setTy( ( muonState.y() - PVPos[1] ) / ( muonState.z() - PVPos[2] ) ); + + // Set the uncertainty on ty to just come from the y uncertainty from the muon stations + auto cov = muonState.covariance(); + cov(3, 3) = std::pow( muonState.ty(), 2 ) * ( cov(1, 1) / ( muonState.y() - PVPos[1] ) ); + muonState.setCovariance( cov ); + // -- Now call the extrapolator and iterate until we have the desired accuracy. double tolerance = 0.5 ; // [mm] Gaudi::TrackMatrix jacobian ; - LHCb::State dummyState = muonState; - for(int i=0; i<10 && std::abs(deltaX) > tolerance ; ++i) { - dummyState = muonState ; - StatusCode sc = m_extrapolator->propagate( dummyState, PVPos[2], &jacobian) ; - if(!sc) return StatusCode::FAILURE; + veloState = muonState; + auto sc = m_extrapolator->propagate( veloState, PVPos[2], &jacobian); + if ( sc.isFailure() ) { + return StatusCode::FAILURE; + } dXdQoP = jacobian(0,4) ; - deltaX = -(dummyState.x() - PVPos[0]) ; + deltaX = -(veloState.x() - PVPos[0]) ; double deltaQoP = deltaX / dXdQoP ; muonState.setQOverP( muonState.qOverP() + deltaQoP ) ; } - finalState = muonState; - return StatusCode::SUCCESS; + veloState.setLocation( LHCb::State::ClosestToBeam ); + muonState.setLocation( LHCb::State::Muon ); + + return m_extrapolator->propagate( veloState, PVPos[2], &jacobian); } diff --git a/Tr/TrackTools/src/MuonSeeding.h b/Tr/TrackTools/src/MuonSeeding.h index 394c545444e..0aafbb802a2 100644 --- a/Tr/TrackTools/src/MuonSeeding.h +++ b/Tr/TrackTools/src/MuonSeeding.h @@ -40,6 +40,7 @@ public: /// Standard constructor using GaudiAlgorithm::GaudiAlgorithm; + MuonSeeding( const std::string& name, ISvcLocator* pSvcLocator ); StatusCode initialize() override; ///< Algorithm initialization StatusCode execute() override; ///< Algorithm execution @@ -47,20 +48,17 @@ private: // -- Methods StatusCode fillPVs(std::vector<double>& PVPos); - void fillMuonStubInfo(LHCb::Track& track, const MuonTrack& muTrack) const; - StatusCode iterateToPV(LHCb::Track* track, LHCb::State& finalState, const std::vector<double>& PVPos, double qOverP); + StatusCode iterateToPV(LHCb::Track* track, LHCb::State& muonState, LHCb::State& veloState, + const std::vector<double>& PVPos, double qOverP); // -- Properties - Gaudi::Property<std::string> m_outputLoc { this, "Output", "Rec/"+name()+"/Tracks" }; - Gaudi::Property<std::string> m_pvLoc { this, "PVLocation", LHCb::RecVertexLocation::Primary }; - Gaudi::Property<std::string> m_extrapolatorName { this, "Extrapolator", "TrackMasterExtrapolator" }; - Gaudi::Property<bool> m_fillMuonStubInfo { this, "FillMuonStubInfo", false }; - Gaudi::Property<std::string> m_trackToolName { this, "ToolName", "MuonNNetRec" }; + Gaudi::Property<bool> m_fitTracks { this, "FitTracks", true }; + Gaudi::Property<std::string> m_outputLoc { this, "Output", "Rec/MuonSeeding/Tracks" }; // -- Tools - IMuonTrackRec* m_trackTool = nullptr; - IMuonTrackMomRec* m_momentumTool = nullptr; - ITrackExtrapolator* m_extrapolator = nullptr; - ITrackFitter* m_trackFitter = nullptr; + ToolHandle<IMuonTrackRec> m_trackTool { "MuonNNetRec/MuonRecTool", this }; + ToolHandle<ITrackFitter> m_trackFitter { "TrackMasterFitter/Fitter", this }; + ToolHandle<IMuonTrackMomRec> m_momentumTool { "MuonTrackMomRec/MuonMomTool", this }; + ToolHandle<ITrackExtrapolator> m_extrapolator { "TrackMasterExtrapolator/Extrapolator", this }; }; #endif // GETMUONTRACK_H -- GitLab From 7e3350d02eb597780acd95e0ac65e992503949c3 Mon Sep 17 00:00:00 2001 From: Olli Lupton <oliver.lupton@cern.ch> Date: Thu, 14 Dec 2017 14:37:47 +0100 Subject: [PATCH 4/5] Revert whitespace change. --- Tr/TrackTools/src/TrackInitFit.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Tr/TrackTools/src/TrackInitFit.h b/Tr/TrackTools/src/TrackInitFit.h index e47cc259e65..dd66d9b3c00 100644 --- a/Tr/TrackTools/src/TrackInitFit.h +++ b/Tr/TrackTools/src/TrackInitFit.h @@ -33,7 +33,7 @@ public: // and fit return m_fitTrack->operator()(track, pid); } - + StatusCode operator() ( std::vector<std::reference_wrapper<LHCb::Track>>& tracks, const LHCb::Tr::PID& pid -- GitLab From 865b224d1d51ca38f33e5f14b09f0124880ec6df Mon Sep 17 00:00:00 2001 From: Ricardo Vazquez Gomez <ricardo.vazquez.gomez@cern.ch> Date: Tue, 20 Feb 2018 09:02:03 +0000 Subject: [PATCH 5/5] change ->fit by ->operator() to be thread safe --- Tr/TrackTools/src/MuonSeeding.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Tr/TrackTools/src/MuonSeeding.cpp b/Tr/TrackTools/src/MuonSeeding.cpp index 4024a48241a..5eac7293983 100644 --- a/Tr/TrackTools/src/MuonSeeding.cpp +++ b/Tr/TrackTools/src/MuonSeeding.cpp @@ -92,7 +92,7 @@ StatusCode MuonSeeding::execute() if ( m_fitTracks ) { // Try and improve the covariance information - auto sc = m_trackFitter->fit( *track ); + auto sc = m_trackFitter->operator()( *track ); if ( sc.isFailure() ) { if ( msgLevel( MSG::WARNING ) ) { warning() << "Track fit failed" << endmsg; -- GitLab