BridgingWireHitRelationFilter.cc

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
#include <tracking/trackFindingCDC/filters/wireHitRelation/BridgingWireHitRelationFilter.h>
 
#include <tracking/trackingUtilities/filters/base/RelationFilter.icc.h>
 
#include <tracking/trackingUtilities/eventdata/hits/CDCWireHit.h>
 
#include <cdc/topology/CDCWireTopology.h>
#include <cdc/topology/CDCWire.h>
 
#include <tracking/trackingUtilities/utilities/Relation.h>
#include <tracking/trackingUtilities/utilities/VectorRange.h>
#include <tracking/trackingUtilities/utilities/Functional.h>
#include <tracking/trackingUtilities/utilities/StringManipulation.h>
 
#include <framework/core/ModuleParamList.templateDetails.h>
 
#include <vector>
#include <string>
#include <cassert>
 
using namespace Belle2;
using namespace CDC;
using namespace TrackFindingCDC;
using namespace TrackingUtilities;
 
template class TrackingUtilities::RelationFilter<const CDCWireHit>;
 
BridgingWireHitRelationFilter::BridgingWireHitRelationFilter() = default;
 
BridgingWireHitRelationFilter::~BridgingWireHitRelationFilter() = default;
 
void BridgingWireHitRelationFilter::exposeParameters(ModuleParamList* moduleParamList,
                                                     const std::string& prefix)
{
  moduleParamList
  ->addParameter(prefixed(prefix, "missingPrimaryNeighborThresholds"),
                 m_param_missingPrimaryNeighborThresholdMap,
                 "Map of o'clock directions to number of missing drift cells "
                 "in the primary neighborhood to trigger the inclusion of secondary neighbors "
                 "in that very o'clock direction",
                 m_param_missingPrimaryNeighborThresholdMap);
}
 
void BridgingWireHitRelationFilter::initialize()
{
  Super::initialize();
  for (short oClockDirection = 0; oClockDirection < 12; oClockDirection++) {
    m_missingPrimaryNeighborThresholds[oClockDirection] = 3;
    if (m_param_missingPrimaryNeighborThresholdMap.count(oClockDirection)) {
      m_missingPrimaryNeighborThresholds[oClockDirection] =
        m_param_missingPrimaryNeighborThresholdMap[oClockDirection];
    }
  }

  for (short oClockDirection : {5, 6, 7, 4, 8, 3, 9, 2, 10, 1, 0, 11}) {
    if (m_missingPrimaryNeighborThresholds[oClockDirection] < 3) {
      m_consideredSecondaryNeighbors.push_back(oClockDirection);
    }
  }
}
 
void BridgingWireHitRelationFilter::prepare(const std::vector<CDCWireHit*>& wireHits)
{
  m_preparedWires.clear();
  m_preparedWires.reserve(wireHits.size());
  for (const CDCWireHit* wireHit : wireHits) {
    m_preparedWires.push_back(&wireHit->getWire());
  }
  m_preparedWireHitsData = wireHits.data();
  m_preparedWireHitsSize = wireHits.size();
}
 
std::vector<CDCWireHit*> BridgingWireHitRelationFilter::getPossibleTos(
  CDCWireHit* from,
  const std::vector<CDCWireHit*>& wireHits) const
{
  assert(std::is_sorted(wireHits.begin(), wireHits.end(), LessOf<Deref>()) &&
         "Expected wire hits to be sorted");
 
  // Use the wires precomputed by prepare() when called with the prepared vector.
  // The wire hits are sorted by the address of their wire such that a search in the
  // contiguous wire array gives the same ranges without dereferencing the wire hits.
  const bool prepared =
    wireHits.data() == m_preparedWireHitsData and wireHits.size() == m_preparedWireHitsSize;
 
  auto findWireHitRange = [&](const CDCWire * neighborWire) -> ConstVectorRange<CDCWireHit*> {
    if (prepared)
    {
      const auto itRange = std::equal_range(m_preparedWires.begin(), m_preparedWires.end(), neighborWire);
      return {
        wireHits.begin() + (itRange.first - m_preparedWires.begin()),
        wireHits.begin() + (itRange.second - m_preparedWires.begin())};
    }
    return ConstVectorRange<CDCWireHit*>{
      std::equal_range(wireHits.begin(), wireHits.end(), neighborWire, LessOf<Deref>())
    };
  };
 
  // Stack buffer for the wire neighbors - at most 6 primary
  // respectively 12 secondary neighbors can be collected
  std::array<std::pair<const CDCWire*, int>, 12> wireNeighbors;
  int nWireNeighbors = 0;
 
  std::vector<CDCWireHit*> wireHitNeighbors;
  wireHitNeighbors.reserve(12);
 
  std::array<short, 12> missingPrimaryNeighbor = {0};
 
  const CDCWireTopology& wireTopology = CDCWireTopology::getInstance();
 
  const CDCWire& wire = from->getWire();
 
  // Analyse primary neighborhood - sorted such that the wire hits relations are most likely sorted.
  for (short oClockDirection : {5, 7, 3, 9, 1, 11}) {
    MayBePtr<const CDCWire> neighborWire = wireTopology.getPrimaryNeighbor(oClockDirection, wire);
    if (neighborWire) wireNeighbors[nWireNeighbors++] = {neighborWire, oClockDirection};
  }
 
  std::sort(wireNeighbors.begin(), wireNeighbors.begin() + nWireNeighbors);
 
  for (std::pair<const CDCWire*, int> wireAndOClockDirection :
       asRange(wireNeighbors.begin(), wireNeighbors.begin() + nWireNeighbors)) {
    const CDCWire* neighborWire = wireAndOClockDirection.first;
    int oClockDirection = wireAndOClockDirection.second;
 
    ConstVectorRange<CDCWireHit*> wireHitRange = findWireHitRange(neighborWire);
    if (wireHitRange.empty()) {
      int ccwOClockDirection = oClockDirection - 1;
      int cwOClockDirection = oClockDirection == 11 ? 0 : oClockDirection + 1;
      ++missingPrimaryNeighbor[ccwOClockDirection];
      ++missingPrimaryNeighbor[oClockDirection];
      ++missingPrimaryNeighbor[cwOClockDirection];
    }
    wireHitNeighbors.insert(wireHitNeighbors.end(), wireHitRange.begin(), wireHitRange.end());
  }
 
  size_t nPrimaryWireHitNeighbors = wireHitNeighbors.size();
  nWireNeighbors = 0;
 
  // Analyse secondary neighborhood
  for (short oClockDirection : m_consideredSecondaryNeighbors) {
    MayBePtr<const CDCWire> neighborWire = wireTopology.getSecondaryNeighbor(oClockDirection, wire);
    if (not neighborWire) continue;
    if (missingPrimaryNeighbor[oClockDirection] <
        m_missingPrimaryNeighborThresholds[oClockDirection])
      continue;
    wireNeighbors[nWireNeighbors++] = {neighborWire, oClockDirection};
  }
 
  std::sort(wireNeighbors.begin(), wireNeighbors.begin() + nWireNeighbors);
 
  for (std::pair<const CDCWire*, int> wireAndOClockDirection :
       asRange(wireNeighbors.begin(), wireNeighbors.begin() + nWireNeighbors)) {
    const CDCWire* neighborWire = wireAndOClockDirection.first;
    ConstVectorRange<CDCWireHit*> wireHitRange = findWireHitRange(neighborWire);
    wireHitNeighbors.insert(wireHitNeighbors.end(), wireHitRange.begin(), wireHitRange.end());
  }

  std::inplace_merge(wireHitNeighbors.begin(),
                     wireHitNeighbors.begin() + nPrimaryWireHitNeighbors,
                     wireHitNeighbors.end(),
                     std::less<CDCWireHit*>());
 
  return wireHitNeighbors;
}