CDCCKFStateCreator.h

/**************************************************************************
 * 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.                  *
 **************************************************************************/
#pragma once
 
#include <tracking/trackingUtilities/findlets/base/Findlet.h>
 
#include <tracking/trackingUtilities/numerics/EForwardBackward.h>
#include <tracking/ckf/general/utilities/SearchDirection.h>
 
#include <tracking/ckf/cdc/entities/CDCCKFState.h>
#include <tracking/ckf/cdc/entities/CDCCKFPath.h>
 
#include <cdc/topology/CDCWire.h>
#include <cdc/topology/CDCWireTopology.h>
 
#include <tracking/trackingUtilities/utilities/StringManipulation.h>
#include <framework/core/ModuleParamList.h>
#include <tracking/trackingUtilities/numerics/Angle.h>
 
#include <Math/Vector3D.h>
 
 
namespace Belle2 {

  class CDCCKFStateCreator
    : public TrackingUtilities::Findlet<CDCCKFState, const CDCCKFState,
      const TrackingUtilities::CDCWireHit* const > {

    using Super = TrackingUtilities::Findlet<CDCCKFState, const CDCCKFState, const TrackingUtilities::CDCWireHit* const>;

    struct CDCCKFWireHitCache {
      int     icLayer;
      double  phi;
    };
 
 
  public:

    void setMaximalLayerJump(int maximalLayerJump) { m_maximalLayerJump = maximalLayerJump; }

    void setMaximalLayerJumpBackwardSeed(int maximalLayerJumpBackwardSeed)
    {
      m_maximalLayerJump_backwardSeed = maximalLayerJumpBackwardSeed;
    }

    void setMaximalDeltaPhi(double maximalDeltaPhi) { m_maximalDeltaPhi = maximalDeltaPhi; }

    void setHitFindingDirection(const std::string& hitFindingDirection)
    {
      m_param_writeOutDirectionAsString = hitFindingDirection;
    }

    void exposeParameters(ModuleParamList* moduleParamList, const std::string& prefix) override
    {
      moduleParamList->addParameter(TrackingUtilities::prefixed(prefix, "maximalLayerJump"),
                                    m_maximalLayerJump, "Maximal jump over N layers", m_maximalLayerJump);
      moduleParamList->addParameter(TrackingUtilities::prefixed(prefix, "maximalLayerJumpBackwardSeed"),
                                    m_maximalLayerJump_backwardSeed, "Maximal jump over N layers", m_maximalLayerJump_backwardSeed);
      moduleParamList->addParameter(TrackingUtilities::prefixed(prefix, "maximalDeltaPhi"),
                                    m_maximalDeltaPhi, "Maximal distance in phi between wires for Z=0 plane", m_maximalDeltaPhi);
      moduleParamList->addParameter(TrackingUtilities::prefixed(prefix, "hitFindingDirection"),
                                    m_param_writeOutDirectionAsString, "Start from innermost/outermost CDC layers", m_param_writeOutDirectionAsString);
 
      moduleParamList->addParameter(TrackingUtilities::prefixed(prefix, "firstActiveCDCLayer"),
                                    m_firstActiveCDCLayer, "First CDC layer which is active (i.e. considered for layer jumps)", m_firstActiveCDCLayer);
    }

    void beginEvent() override
    {
      Super::beginEvent();
      m_wireHitCache.clear();
 
      // Determine direction of track building
      m_param_writeOutDirection = fromString(m_param_writeOutDirectionAsString);
 
      if (m_param_writeOutDirection == TrackingUtilities::EForwardBackward::c_Forward) {
        doForward = true;
      } else if (m_param_writeOutDirection == TrackingUtilities::EForwardBackward::c_Backward) {
        doForward = false;
      } else {
        B2FATAL("CDCCKFStateCreator: No valid direction specified. Please use forward/backward.");
      }
    }

    void apply(std::vector<CDCCKFState>& nextStates, const CDCCKFPath& path,
               const std::vector<const TrackingUtilities::CDCWireHit*>& wireHits) override
    {
      // TODO: as we do not need any information on the current state (track state) of the path, we could in principle
      // TODO: precalculate everything in here
 
      // Create cache over wirehits, if empty:
      if (m_wireHitCache.empty()) {
        const size_t nHits = wireHits.size();
        m_wireHitCache.reserve(nHits);
        for (auto  hitPtr : wireHits) {
          // to speed things up, don't consider background/taken hits at all (and not just in the loop below).
          // I can't just remove them from the list, otherwise the relation to the wireHits is broken
          // so set the layer index to a high number.
          if (hitPtr->getAutomatonCell().hasBackgroundFlag() || hitPtr->getAutomatonCell().hasTakenFlag()) {
            m_wireHitCache.push_back(CDCCKFWireHitCache{99999, 0.});
          } else {
            m_wireHitCache.push_back(CDCCKFWireHitCache{hitPtr->getWire().getICLayer(), hitPtr->getRefPos2D().Phi()});
          }
        }
      }
 
      // Cache last-on-the-path state info too:
      const auto& lastState = path.back();
      double lastPhi = 0;
      double lastICLayer = -1;
      if (lastState.isSeed()) {
        if (doForward) {
          lastICLayer = m_firstActiveCDCLayer; // if SL0 and SL1 turned off this 14, nominal 0
        } else {
          const auto& wireTopology = CDC::CDCWireTopology::getInstance();
          const auto& wires = wireTopology.getWires();
          const float maxForwardZ = wires.back().getForwardZ();     // 157.615
          const float maxBackwardZ = wires.back().getBackwardZ();   // -72.0916
 
          const ROOT::Math::XYZVector seedPos(lastState.getSeed()->getPositionSeed());
          const float seedPosZ = seedPos.z();
 
          if (seedPosZ < maxForwardZ && seedPosZ > maxBackwardZ) {
            lastICLayer = 56;
          } else {
            // do straight extrapolation of seed momentum to CDC outer walls
            ROOT::Math::XYZVector seedMomZOne(lastState.getSeed()->getMomentumSeed());
            seedMomZOne = seedMomZOne / seedMomZOne.z();
            // const float maxZ = seedPosZ > 0 ? maxForwardZ : maxBackwardZ;
            // const ROOT::Math::XYZVector extrapolatedPos = seedPos - seedMom / seedMom.R() * (seedPosZ - maxZ);
 
            // find closest iCLayer
            float minDist = 99999;
            for (const auto& wire : wires) {
              const float maxZ = seedPosZ > 0 ? wire.getForwardZ() : wire.getBackwardZ();
              const ROOT::Math::XYZVector extrapolatedPos = seedPos - seedMomZOne * (seedPosZ - maxZ);
 
              const auto distance = wire.getDistance(extrapolatedPos);
              if (distance < minDist) {
                minDist = distance;
                lastICLayer = wire.getICLayer();
              }
            }
            B2DEBUG(29, lastICLayer << " (d=" << minDist << ")");
          }
        }
      } else {
        lastPhi = lastState.getWireHit()->getRefPos2D().Phi();
        lastICLayer = lastState.getWireHit()->getWire().getICLayer();
      }
 
      // Get sorted vector of wireHits on the path for faster search
      std::vector<const TrackingUtilities::CDCWireHit*> wireHitsOnPath;
      for (auto const& state : path) {
        if (! state.isSeed()) {
          wireHitsOnPath.push_back(state.getWireHit());
        }
      }
      std::sort(wireHitsOnPath.begin(), wireHitsOnPath.end());
 
      size_t nHits = wireHits.size();
      for (size_t i = 0; i < nHits; i++) {
        // adjust direction of loop (minimal speed gain)
        int idx = doForward ? i : nHits - i - 1;
 
        const auto iCLayer =  m_wireHitCache[idx].icLayer; // wireHit->getWire().getICLayer();
        if (m_param_writeOutDirection == TrackingUtilities::EForwardBackward::c_Backward && lastState.isSeed()) {
          if (std::abs(lastICLayer - iCLayer) > m_maximalLayerJump_backwardSeed) {
            continue;
          }
        } else if (std::abs(lastICLayer - iCLayer) > m_maximalLayerJump) {
          continue;
        }
 
        if (! lastState.isSeed()) {
          double deltaPhi = TrackingUtilities::AngleUtil::normalised(lastPhi - m_wireHitCache[idx].phi);
          if (fabs(deltaPhi)  > m_maximalDeltaPhi)  {
            continue;
          }
        }
 
        const TrackingUtilities::CDCWireHit* wireHit = wireHits[idx];
 
        if (std::binary_search(wireHitsOnPath.begin(), wireHitsOnPath.end(), wireHit)) {
          continue;
        }
 
        nextStates.emplace_back(wireHit);
      }
    }
 
  private:
    int m_maximalLayerJump = 2;
    int m_firstActiveCDCLayer = 0;
    int m_maximalLayerJump_backwardSeed = 3;
    double m_maximalDeltaPhi =  TMath::Pi() / 8;
    std::string m_param_writeOutDirectionAsString = "forward";
    TrackingUtilities::EForwardBackward m_param_writeOutDirection = TrackingUtilities::EForwardBackward::c_Unknown;
    bool doForward = true;

    std::vector<CDCCKFWireHitCache> m_wireHitCache = {};
 
  };

}