REC: Silicon tracking

Reconstruction/CMakeLists.txt

@@ -4,3 +4,4 @@ add_subdirectory(PFA)
 add_subdirectory(SiliconTracking)
 add_subdirectory(Tracking)
 add_subdirectory(RecGenfitAlg)
+add_subdirectory(RecAssociationMaker)

Reconstruction/RecAssociationMaker/CMakeLists.txt

@@ -0,0 +1,20 @@
+# Modules
+gaudi_add_module(RecAssociationMaker
+                 SOURCES src/TrackParticleRelationAlg.cpp
+
+                 LINK GearSvc
+                      EventSeeder
+                      TrackSystemSvcLib
+                      DataHelperLib
+                      KiTrackLib
+                      Gaudi::GaudiKernel
+                      k4FWCore::k4FWCore
+                      ${GEAR_LIBRARIES}
+                      ${GSL_LIBRARIES}
+                      ${LCIO_LIBRARIES}
+)
+install(TARGETS RecAssociationMaker
+  EXPORT CEPCSWTargets
+  RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
+  LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib
+  COMPONENT dev)

Reconstruction/RecAssociationMaker/src/TrackParticleRelationAlg.cpp

@@ -0,0 +1,150 @@
+#include "TrackParticleRelationAlg.h"
+
+DECLARE_COMPONENT(TrackParticleRelationAlg)
+
+TrackParticleRelationAlg::TrackParticleRelationAlg(const std::string& name, ISvcLocator* svcLoc)
+: GaudiAlgorithm(name, svcLoc){
+  declareProperty("MCParticleCollection", m_inMCParticleColHdl, "Handle of the Input MCParticle collection");
+}
+
+StatusCode TrackParticleRelationAlg::initialize() {
+  for(auto name : m_inTrackCollectionNames) {
+    m_inTrackColHdls.push_back(new DataHandle<edm4hep::TrackCollection> (name, Gaudi::DataHandle::Reader, this));
+    m_outColHdls.push_back(new DataHandle<edm4hep::MCRecoTrackParticleAssociationCollection> (name+"ParticleAssociation", Gaudi::DataHandle::Writer, this));
+  }
+
+  for(unsigned i=0; i<m_inAssociationCollectionNames.size(); i++) {
+    m_inAssociationColHdls.push_back(new DataHandle<edm4hep::MCRecoTrackerAssociationCollection> (m_inAssociationCollectionNames[i], Gaudi::DataHandle::Reader, this));
+  }
+
+  if(m_inAssociationColHdls.size()==0) {
+    warning() << "no Association Collection input" << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  m_nEvt = 0;
+  return GaudiAlgorithm::initialize();
+}
+
+StatusCode TrackParticleRelationAlg::execute() {
+  info() << "start Event " << m_nEvt << endmsg;
+
+  // MCParticle
+  const edm4hep::MCParticleCollection* mcpCol = nullptr;
+  try {
+    mcpCol = m_inMCParticleColHdl.get();
+    for (auto mcp : *mcpCol) {
+      debug() << "id=" << mcp.id() << " pdgid=" << mcp.getPDG() << " charge=" << mcp.getCharge() << " genstat=" << mcp.getGeneratorStatus() << " simstat=" << mcp.getSimulatorStatus()
+	      << " p=" << mcp.getMomentum() << endmsg;
+      int nparents = mcp.parents_size();
+      int ndaughters = mcp.daughters_size();
+      for (int i=0; i<nparents; i++) {
+	debug() << "<<<<<<" << mcp.getParents(i).id() << endmsg;
+      }
+      for (int i=0; i<ndaughters; i++) {
+        debug() << "<<<<<<" << mcp.getDaughters(i).id() << endmsg;
+      }
+    }
+  }
+  catch ( GaudiException &e ) {
+    debug() << "Collection " << m_inMCParticleColHdl.fullKey() << " is unavailable in event " << m_nEvt << endmsg;
+  }
+  if(!mcpCol) {
+    warning() << "pass this event because MCParticle collection can not read" << endmsg; 
+    return StatusCode::SUCCESS;
+  }
+
+  // Prepare map from hit to MCParticle
+  std::map<edm4hep::TrackerHit, edm4hep::MCParticle> mapHitParticle;
+  debug() << "reading Association" << endmsg;
+  for (auto hdl : m_inAssociationColHdls) {
+    const edm4hep::MCRecoTrackerAssociationCollection* assCol = nullptr;
+    try {
+      assCol = hdl->get();
+    }
+    catch ( GaudiException &e ) {
+      debug() << "Collection " << hdl->fullKey() << " is unavailable in event " << m_nEvt << endmsg;
+      return StatusCode::FAILURE;
+    }
+    for (auto ass: *assCol) {
+      auto hit = ass.getRec();
+      auto particle = ass.getSim().getMCParticle();
+      mapHitParticle[hit] = particle;
+    }
+  }
+
+  // Handle all input TrackCollection
+  for (unsigned icol=0; icol<m_inTrackColHdls.size(); icol++) {
+    auto hdl = m_inTrackColHdls[icol];
+
+    const edm4hep::TrackCollection* trkCol = nullptr;
+    try {
+      trkCol = hdl->get();
+    }
+    catch ( GaudiException &e ) {
+      debug() << "Collection " << hdl->fullKey() << " is unavailable in event " << m_nEvt << endmsg;
+    }
+
+    auto outCol = m_outColHdls[icol]->createAndPut();
+    if(!outCol) {
+      error() << "Collection " << m_outColHdls[icol]->fullKey() << " cannot be created" << endmsg;
+      return StatusCode::FAILURE;
+    }
+
+    if(trkCol) {
+      std::map<edm4hep::MCParticle, std::vector<edm4hep::TrackerHit> > mapParticleHits;
+
+      for (auto track: *trkCol) {
+	std::map<edm4hep::MCParticle, int> mapParticleNHits;
+
+	// Count hits related to MCParticle
+	int nhits = track.trackerHits_size();
+	for (int ihit=0; ihit<nhits; ihit++) {
+	  auto hit = track.getTrackerHits(ihit);
+	  auto itHP  = mapHitParticle.find(hit);
+	  if (itHP==mapHitParticle.end()) {
+	    warning() << "track " << track.id() << "'s hit " << hit.id() << "not in association list, will be ignored" << endmsg;
+	  }
+	  else {
+	    auto particle = itHP->second;
+	    if(!particle.isAvailable()) continue;
+	    debug() << "track " << track.id() << "'s hit " << hit.id() << " link to MCParticle " << particle.id() << endmsg;
+	    auto itPN = mapParticleNHits.find(particle);
+	    if (itPN!=mapParticleNHits.end()) itPN->second++;
+	    else                             mapParticleNHits[particle] = 1;
+
+	    if (msgLevel(MSG::DEBUG)) mapParticleHits[particle].push_back(hit);
+	  }
+	}
+	debug() << "Total " << mapParticleNHits.size() << " particles link to the track " << track.id() << endmsg;
+
+	// Save to collection
+	for (std::map<edm4hep::MCParticle, int>::iterator it=mapParticleNHits.begin(); it!=mapParticleNHits.end(); it++) {
+	  auto ass = outCol->create();
+	  ass.setWeight(it->second);
+	  ass.setRec(track);
+	  ass.setSim(it->first);
+	}
+      }
+
+      if (msgLevel(MSG::DEBUG)) {
+	for (std::map<edm4hep::MCParticle, std::vector<edm4hep::TrackerHit> >::iterator it=mapParticleHits.begin(); it!=mapParticleHits.end(); it++) {
+	  auto particle = it->first;
+	  auto hits     = it->second;
+	  debug() << "=== MCPaticle ===" << particle << endmsg;
+	  for (auto hit : hits) {
+	    debug() << hit << endmsg;
+	  }
+	}
+      }
+    }
+  }
+
+  m_nEvt++;
+  return StatusCode::SUCCESS;
+}
+
+StatusCode TrackParticleRelationAlg::finalize() {
+
+  return StatusCode::SUCCESS;
+}

Reconstruction/RecAssociationMaker/src/TrackParticleRelationAlg.h

@@ -0,0 +1,38 @@
+#ifndef TrackParticleRelationAlg_h
+#define TrackParticleRelationAlg_h 1
+
+#include "k4FWCore/DataHandle.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/TrackCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/MCRecoTrackParticleAssociationCollection.h"
+
+class TrackParticleRelationAlg : public GaudiAlgorithm {
+ public:
+  // Constructor of this form must be provided
+  TrackParticleRelationAlg( const std::string& name, ISvcLocator* pSvcLocator );
+
+  // Three mandatory member functions of any algorithm
+  StatusCode initialize() override;
+  StatusCode execute() override;
+  StatusCode finalize() override;
+
+ private:
+  // input MCParticle
+  DataHandle<edm4hep::MCParticleCollection>                              m_inMCParticleColHdl{"MCParticle", Gaudi::DataHandle::Reader, this};
+  // input Tracks to make relation
+  std::vector<DataHandle<edm4hep::TrackCollection>* >                    m_inTrackColHdls;
+  Gaudi::Property<std::vector<std::string> >                             m_inTrackCollectionNames{this, "TrackList", {"SiTracks"}};
+  // input TrackerAssociation to link TrackerHit and SimTrackerHit
+  std::vector<DataHandle<edm4hep::MCRecoTrackerAssociationCollection>* > m_inAssociationColHdls;
+  Gaudi::Property<std::vector<std::string> >                             m_inAssociationCollectionNames{this, "TrackerAssociationList", {"VXDTrackerHitAssociation",
+        "SITTrackerHitAssociation", "SETTrackerHitAssociation", "FTDTrackerHitAssociation"}};
+
+  // output TrackParticleAssociation
+  std::vector<DataHandle<edm4hep::MCRecoTrackParticleAssociationCollection>* > m_outColHdls;
+
+  int m_nEvt;
+};
+#endif

Reconstruction/SiliconTracking/CMakeLists.txt

@@ -7,6 +7,7 @@ gaudi_add_module(SiliconTracking
                          src/TrackSubsetAlg.cpp
                  LINK GearSvc
                       EventSeeder
+                      TrackingLib
                       TrackSystemSvcLib
                       DataHelperLib
                       KiTrackLib

Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp

@@ -4,7 +4,7 @@
 #include "DataHelper/Navigation.h"
 
 #include "edm4hep/TrackerHit.h"
-#include "edm4hep/TrackerHit.h"
+#include "edm4hep/TrackerHitPlane.h"
 #include "edm4hep/Track.h"
 #if __has_include("edm4hep/EDM4hepVersion.h")
 #include "edm4hep/EDM4hepVersion.h"
@@ -74,12 +74,16 @@ StatusCode ForwardTrackingAlg::initialize(){
   // can be printed. As this is mainly used for debugging it is not a steerable parameter.
   //if( _useCED )MarlinCED::init(this) ;    //CED
 
+  // Set up the track fit tool
+  m_fitTool = ToolHandle<ITrackFitterTool>(m_fitToolName.value());
+
   if(m_dumpTime){
     NTuplePtr nt1(ntupleSvc(), "MyTuples/Time"+name());
     if ( !nt1 ) {
       m_tuple = ntupleSvc()->book("MyTuples/Time"+name(),CLID_ColumnWiseTuple,"Tracking time");
       if ( 0 != m_tuple ) {
 	m_tuple->addItem ("timeTotal",  m_timeTotal ).ignore();
+	m_tuple->addItem ("timeKalman", m_timeKalman ).ignore();
       }
       else {
 	fatal() << "Cannot book MyTuples/Time"+name() <<endmsg;
@@ -306,7 +310,9 @@ StatusCode ForwardTrackingAlg::execute(){
       _map_sector_hits[ ftdHit->getSector() ].push_back( ftdHit );         
     }
   }
-
+
+  if(m_dumpTime&&m_tuple) m_timeKalman = 0;
+
   if( !_map_sector_hits.empty() ){
     /**********************************************************************************************/
     /*                Check if no sector is overflowing with hits                                 */
@@ -661,7 +667,7 @@ StatusCode ForwardTrackingAlg::execute(){
     debug() << "\t\t---Save Tracks---" << endmsg ;
 
     //auto trkCol = _outColHdl.createAndPut();
-    
+
     for (unsigned int i=0; i < tracks.size(); i++){
       FTDTrack* myTrack = dynamic_cast< FTDTrack* >( tracks[i] );
 
@@ -942,7 +948,7 @@ bool ForwardTrackingAlg::setCriteria( unsigned round ){
 }
 
 void ForwardTrackingAlg::finaliseTrack( edm4hep::MutableTrack* trackImpl ){
-     
+  auto stopwatch = TStopwatch();
   Fitter fitter( trackImpl , _trkSystem );
 
   //trackImpl->trackStates().clear();
@@ -1024,6 +1030,9 @@ void ForwardTrackingAlg::finaliseTrack( edm4hep::MutableTrack* trackImpl ){
   trackImpl->addToSubDetectorHitNumbers(hitNumbers[UTIL::ILDDetID::SET]);
   trackImpl->addToSubDetectorHitNumbers(hitNumbers[UTIL::ILDDetID::ETD]);
 #endif
+
+  if(m_dumpTime&&m_tuple) m_timeKalman += stopwatch.RealTime()*1000;
+
   return;
 }
 

Reconstruction/SiliconTracking/src/ForwardTrackingAlg.h

@@ -14,6 +14,7 @@
 
 #include "edm4hep/TrackerHit.h"
 #include "TrackSystemSvc/IMarlinTrkSystem.h"
+#include "Tracking/ITrackFitterTool.h"
 //#include "gear/BField.h"
 
 #include "KiTrack/Segment.h"
@@ -239,6 +240,7 @@ class ForwardTrackingAlg : public GaudiAlgorithm {
 
   NTuple::Tuple*       m_tuple;
   NTuple::Item<float>  m_timeTotal;
+  NTuple::Item<float>  m_timeKalman;
 
   /** B field in z direction */
   double _Bz;
@@ -261,6 +263,7 @@ class ForwardTrackingAlg : public GaudiAlgorithm {
   Gaudi::Property<std::vector<float> > _critMinimaInit{this, "CriteriaMin", {} };
   Gaudi::Property<std::vector<float> > _critMaximaInit{this, "CriteriaMax", {} };
   Gaudi::Property<bool>   m_dumpTime{this, "DumpTime", true};
+  Gaudi::Property<std::string> m_fitToolName{this, "FitterTool", "KalTestTool/KalTest110"};
 
   std::map<std::string, std::vector<float> > _critMinima;
   std::map<std::string, std::vector<float> > _critMaxima;
@@ -291,6 +294,7 @@ class ForwardTrackingAlg : public GaudiAlgorithm {
   unsigned _nTrackCandidatesPlus;
 
   MarlinTrk::IMarlinTrkSystem* _trkSystem;
+  ToolHandle<ITrackFitterTool> m_fitTool;
 
   /** The quality of the output track collection */
   int _output_track_col_quality ;