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