Newer
Older
TB_Chris / TbAlgorithms / src / .svn / text-base / TbClusterAssociator.cpp.svn-base
  1. // Gaudi
  2. #include "GaudiKernel/PhysicalConstants.h"
  3.  
  4. // Tb/TbEvent
  5. #include "Event/TbTrack.h"
  6. #include "Event/TbCluster.h"
  7.  
  8. // Tb/TbKernel
  9. #include "TbKernel/TbConstants.h"
  10.  
  11. // Local
  12. #include "TbClusterAssociator.h"
  13.  
  14. DECLARE_ALGORITHM_FACTORY(TbClusterAssociator)
  15.  
  16. //=============================================================================
  17. // Standard constructor
  18. //=============================================================================
  19. TbClusterAssociator::TbClusterAssociator(const std::string& name,
  20. ISvcLocator* pSvcLocator)
  21. : TbAlgorithm(name, pSvcLocator) {
  22.  
  23. declareProperty("TrackLocation",
  24. m_trackLocation = LHCb::TbTrackLocation::Default);
  25. declareProperty("ClusterLocation",
  26. m_clusterLocation = LHCb::TbClusterLocation::Default);
  27.  
  28. declareProperty("DUTs", m_duts);
  29.  
  30. declareProperty("UseHits", m_useHits = true);
  31. declareProperty("ReuseClusters", m_reuseClusters = true);
  32. declareProperty("TimeWindow", m_twindow = 200. * Gaudi::Units::ns);
  33. declareProperty("XWindow", m_xwindow = 1. * Gaudi::Units::mm);
  34. declareProperty("MaxChi2", m_maxChi2 = 100.);
  35. }
  36.  
  37. //=============================================================================
  38. // Initialization
  39. //=============================================================================
  40. StatusCode TbClusterAssociator::initialize() {
  41.  
  42. // Initialise the base class.
  43. StatusCode sc = TbAlgorithm::initialize();
  44. if (sc.isFailure()) return sc;
  45. return StatusCode::SUCCESS;
  46. }
  47.  
  48. //=============================================================================
  49. // Main execution
  50. //=============================================================================
  51. StatusCode TbClusterAssociator::execute() {
  52.  
  53. // Grab the tracks.
  54. const LHCb::TbTracks* tracks = getIfExists<LHCb::TbTracks>(m_trackLocation);
  55. if (!tracks) {
  56. return Error("No tracks in " + m_trackLocation);
  57. }
  58.  
  59. for (const auto dut : m_duts) {
  60. // Get the clusters for this plane.
  61. const std::string clusterLocation = m_clusterLocation + std::to_string(dut);
  62. const LHCb::TbClusters* clusters =
  63. getIfExists<LHCb::TbClusters>(clusterLocation);
  64. if (!clusters) continue;
  65. // Loop over the tracks.
  66. for (LHCb::TbTrack* track : *tracks) {
  67. // Skip low quality tracks.
  68. if (track->chi2PerNdof() > m_maxChi2) continue;
  69. const Gaudi::XYZPoint pGlobal = geomSvc()->intercept(track, dut);
  70. const auto pLocal = geomSvc()->globalToLocal(pGlobal, dut);
  71. const double tMin = track->htime() - m_twindow;
  72. const double tMax = track->htime() + m_twindow;
  73. for (LHCb::TbCluster* cluster : *clusters) {
  74. // Assume that the clusters are sorted by time.
  75. if (cluster->htime() < tMin) continue;
  76. if (cluster->htime() > tMax) break;
  77. // Skip used clusters (if requested).
  78. if (!m_reuseClusters && cluster->associated()) continue;
  79. if (m_useHits) {
  80. // Compare the coordinates of the pixels to the track intercept.
  81. if (!match(cluster, pLocal.x(), pLocal.y())) continue;
  82. } else {
  83. // Compare the cluster coordinates to the track intercept.
  84. const double dx = cluster->xloc() - pLocal.x();
  85. const double dy = cluster->yloc() - pLocal.y();
  86. if (fabs(dx) > m_xwindow || fabs(dy) > m_xwindow) continue;
  87. }
  88. // Associate the cluster to the track and tag it as used.
  89. track->addToAssociatedClusters(cluster);
  90. cluster->setAssociated(true);
  91. }
  92. }
  93. }
  94. return StatusCode::SUCCESS;
  95. }
  96.  
  97. //=============================================================================
  98. // Check if the cluster has hits within the search window.
  99. //=============================================================================
  100. bool TbClusterAssociator::match(const LHCb::TbCluster* cluster, const double x,
  101. const double y) const {
  102.  
  103. const unsigned int plane = cluster->plane();
  104. // Loop over hits in the cluster.
  105. for (auto hit : cluster->hits()) {
  106. double xLocal = 0.;
  107. double yLocal = 0.;
  108. geomSvc()->pixelToPoint(hit->scol(), hit->row(), plane, xLocal, yLocal);
  109. const double dx = xLocal - x;
  110. const double dy = yLocal - y;
  111. if (fabs(dx) < m_xwindow && fabs(dy) < m_xwindow) return true;
  112. }
  113. return false;
  114. }