diff --git a/Muon/MuonTrackRec/src/component/MuonCombRec.cpp b/Muon/MuonTrackRec/src/component/MuonCombRec.cpp index 59d5821f3f7c6ef06ca01385e5b48ed714f5538b..2ff18db7bc6b4ed94329791779298ef476ec5f9b 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 5600f80aa7051c974354e86eba702f4a8f525959..0d44f766c4d9ff67c306cb5c3d96cf2d50c35a7b 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 9ad18771522a7b828b08a55fc92f1bd924c7855a..9b6d93fb6aeb4349830b9de1ac20f7b57d66bae1 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 88f0923435a36b00151133ad5ad2825cf77b6aa7..a803cd3b96a8ca60927541e6a969d2bb5d7bd887 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 83e78cf0e70202633a310ab4eaaf05532d09cc5f..08cb6e45932fe9af231ae6fb77ef6643195d8f64 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/Pr/PrMCTools/src/PrMuonStubAssociator.cpp b/Pr/PrMCTools/src/PrMuonStubAssociator.cpp new file mode 100644 index 0000000000000000000000000000000000000000..59ba6f8dbe4ce6dd04c943c32c8cdbf8a41a38d6 --- /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 0000000000000000000000000000000000000000..8ac429d6a7387866c5b865b1aeb4b63763b6ae95 --- /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 new file mode 100644 index 0000000000000000000000000000000000000000..5eac7293983fa75cbca20a1672143293ebebb539 --- /dev/null +++ b/Tr/TrackTools/src/MuonSeeding.cpp @@ -0,0 +1,184 @@ +// Include files + +// local +#include "MuonSeeding.h" +#include "Event/VPCluster.h" + +//----------------------------------------------------------------------------- +// Implementation file for class : MuonSeeding +// +// 2010-09-14 : Michel De Cian +//----------------------------------------------------------------------------- + +// 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 +//============================================================================= +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.retrieve(); + m_extrapolator.retrieve(); + m_trackFitter.retrieve(); + m_momentumTool.retrieve(); + + 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 }; // 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 ); + } + + 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 ( const auto& muonTrack : muonTracks ) { + // -- Get the momentum of the muon track, make a new track + 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->operator()( *track ); + if ( sc.isFailure() ) { + if ( msgLevel( MSG::WARNING ) ) { + warning() << "Track fit failed" << endmsg; + } + + continue; + } + } + + // -- 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; + } + + // -- Set Pattern Reco status and track type, finally fit the track + track->setPatRecStatus( LHCb::Track::PatRecStatus::PatRecIDs ); + track->setType( LHCb::Track::Types::Muon ); + + if ( !m_fitTracks ) { + track->clearStates( ); // remove the state recMomentum created + } + + track->addToStates( veloState ); + track->addToStates( muonState ); + // --------------------------------------------------------------- + + // -- Finally, insert the track into the container! + tracks->insert( track.release() ); + } + + // --------------------------------------------------------------- + + if(msgLevel(MSG::DEBUG)) { + debug() << "Filling " << tracks->size() << " tracks in " << m_outputLoc << endmsg; + } + + // Control flow can stop if we didn't make anything + setFilterPassed( !tracks->empty() ); + + // Put our output on the TES + put( tracks.release(), m_outputLoc ); + + return StatusCode::SUCCESS; +} + +//============================================================================= +// Change the q/p till the track points to the PV (stolen from Wouter) +//============================================================================= +StatusCode MuonSeeding::iterateToPV(LHCb::Track* track, LHCb::State& muonState, LHCb::State& veloState, + const std::vector<double> &PVPos, double qOverP) +{ + 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 ; + + for(int i=0; i<10 && std::abs(deltaX) > tolerance ; ++i) { + veloState = muonState; + auto sc = m_extrapolator->propagate( veloState, PVPos[2], &jacobian); + if ( sc.isFailure() ) { + return StatusCode::FAILURE; + } + dXdQoP = jacobian(0,4) ; + deltaX = -(veloState.x() - PVPos[0]) ; + double deltaQoP = deltaX / dXdQoP ; + muonState.setQOverP( muonState.qOverP() + deltaQoP ) ; + } + + 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 new file mode 100644 index 0000000000000000000000000000000000000000..0aafbb802a20d20c8acfb7b2141016a8d0051e4e --- /dev/null +++ b/Tr/TrackTools/src/MuonSeeding.h @@ -0,0 +1,64 @@ +#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 "Event/State.h" +#include "Event/RecVertex.h" + +#include "TrackInterfaces/ITrackFitter.h" +#include "TrackInterfaces/ITrackExtrapolator.h" + +#include "Kernel/LHCbID.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. + * - 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; + + MuonSeeding( const std::string& name, ISvcLocator* pSvcLocator ); + StatusCode initialize() override; ///< Algorithm initialization + StatusCode execute() override; ///< Algorithm execution + +private: + + // -- Methods + StatusCode fillPVs(std::vector<double>& PVPos); + StatusCode iterateToPV(LHCb::Track* track, LHCb::State& muonState, LHCb::State& veloState, + const std::vector<double>& PVPos, double qOverP); + + // -- Properties + Gaudi::Property<bool> m_fitTracks { this, "FitTracks", true }; + Gaudi::Property<std::string> m_outputLoc { this, "Output", "Rec/MuonSeeding/Tracks" }; + + // -- Tools + 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