NDFinder.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 <string>
#include <vector>
#include <array>
#include <utility>
#include <Math/Vector3D.h>
#include "trg/cdc/NDFinderDefs.h"
#include "trg/cdc/Clusterizend.h"
 
namespace Belle2 {

  struct NDFinderParameters {
    unsigned short minSuperAxial;
    unsigned short minSuperStereo;
    unsigned short iterations;
    unsigned short omegaTrim;
    unsigned short phiTrim;
    bool storeAdditionalReadout;
    std::string axialFile;
    std::string stereoFile;
  };

  struct HitInfo {
    unsigned short hitID;
    unsigned short hitSLID;
    unsigned short hitPrioPos;
    short hitPrioTime;
  };

  class NDFinderTrack {
  public:
    NDFinderTrack(std::array<double, 3> estimatedParameters,
                  SimpleCluster&& cluster,
                  std::vector<ROOT::Math::XYZVector>&& readoutHoughSpace,
                  std::vector<ROOT::Math::XYZVector>&& readoutCluster)
      : m_cluster(std::move(cluster)),
        m_houghSpace(std::move(readoutHoughSpace)),
        m_readoutCluster(std::move(readoutCluster))
    {
      m_omega = estimatedParameters[0];
      m_phi = estimatedParameters[1];
      m_cot = estimatedParameters[2];
    }

    double getOmega() const { return m_omega; }
    double getPhi0() const { return m_phi; }
    double getCot() const { return m_cot; }
    std::vector<unsigned short> getRelatedHits() const { return m_cluster.getClusterHits(); }
    std::vector<ROOT::Math::XYZVector> getHoughSpace() const { return m_houghSpace; }
    std::vector<ROOT::Math::XYZVector> getClusterReadout() const { return m_readoutCluster; }
 
  private:
    double m_omega;
    double m_phi;
    double m_cot;
    SimpleCluster m_cluster;
    std::vector<ROOT::Math::XYZVector> m_houghSpace;
    std::vector<ROOT::Math::XYZVector> m_readoutCluster;
  };
 

  class NDFinder {
  public:
    struct SectorBinning {
      unsigned short omega;
      unsigned short phi;
      unsigned short cot;
      unsigned short nHitIDs;
      unsigned short nPriorityWires;
    };

    struct WireInfo {
      unsigned short relativeWireID;
      unsigned short phiSectorStart;
      unsigned short priorityWire;
    };

    struct ContributionInfo {
      unsigned short hitIndex;
      unsigned short contribution;
      unsigned short superLayer;
      short priorityTime;
    };
 
    // Default constructor
    NDFinder() = default;
 
    // Destructor
    virtual ~NDFinder()
    {
      delete m_hitToSectorIDs;
      delete m_compAxialHitReps;
      delete m_compStereoHitReps;
      delete m_expAxialHitReps;
      delete m_expStereoHitReps;
      delete m_houghSpace;
    }

    void init(const NDFinderParameters& ndFinderParameters);
    void reset();
    void addHit(const HitInfo& hitInfo);
    void findTracks();
    std::vector<NDFinderTrack>* getFinderTracks() { return &m_ndFinderTracks; }
 
    // NDFinder: Internal functions for track finding
  protected:
    void initLookUpArrays();
    void initHitToSectorMap();
    void loadCompressedHitReps(const std::string& fileName, const SectorBinning& compBins, c5array& compHitsToWeights);
    void fillExpandedHitReps(const SectorBinning& compBins, const c5array& compHitsToWeights, c5array& expHitsToWeights);
    void processHitForHoughSpace(const unsigned short hitIdx);
    void writeHitToHoughSpace(const WireInfo& hitInfo, const c5array& expHitsToWeights);
    void runTrackFinding();
    std::vector<SimpleCluster> relateHitsToClusters(std::vector<SimpleCluster>& clusters);
    std::vector<std::vector<unsigned short>> getHitsVsClustersTable(const std::vector<SimpleCluster>& clusters);
    unsigned short getHitContribution(const cell_index& peakCell, const unsigned short hitIdx);
    std::vector<ContributionInfo> extractContributionInfos(const std::vector<unsigned short>& clusterHits);
    int getMaximumHitInSuperLayer(const std::vector<ContributionInfo>& contributionInfos, unsigned short superLayer);
    bool checkHitSuperLayers(const SimpleCluster& cluster);
    std::array<double, 3> calculateCenterOfGravity(const SimpleCluster& cluster);
    std::array<double, 3> getTrackParameterEstimate(const std::array<double, 3>& centerOfGravity);
    std::array<double, 3> transformTrackParameters(const std::array<double, 3>& estimatedParameters);
    static inline double getTrackRadius(double transverseMomentum) { return transverseMomentum * 1e11 / (3e8 * 1.5); }
 
    // NDFinder: Member data stores
  private:
    std::vector<NDFinderTrack> m_ndFinderTracks;
    std::vector<unsigned short> m_hitIDs;
    std::vector<unsigned short> m_hitSLIDs;
    std::vector<unsigned short> m_priorityWirePos;
    std::vector<short> m_priorityWireTime;
    unsigned short m_nHits{0};
    NDFinderParameters m_ndFinderParams;
    Clusterizend m_clusterer;
 
    /*
      Since the CDC wire pattern is repeated 32 times, the hit IDs are stored for 1/32 of the CDC only.
      The total number of 2336 TS corresponds to (41 axial + 32 stereo) * 32.
      The number of track bins is (for example): (omega, phi, cot) = (40, 384, 9)
      Note: The omega dimension here represents just sign(q)/(p_T[GeV/c]) (0.25 -> inf -> -inf -> -0.25 for 0 -> 19.5 -> 39)
    */
 
    // Track segments
    static constexpr unsigned short m_nTS = 2336; 
    static constexpr unsigned short m_nSL = 9; 
    static constexpr unsigned short m_nAxial = 41; 
    static constexpr unsigned short m_nStereo = 32; 
    static constexpr unsigned short m_nPrio = 3; 
 
    // Full Hough space bins
    static constexpr unsigned short m_nOmega = 40; 
    static constexpr unsigned short m_nPhi = 384; 
    static constexpr unsigned short m_nCot = 9; 
 
    // CDC symmetry in phi
    static constexpr unsigned short m_phiGeo = 32; 
 
    // Phi sectors in the CDC
    static constexpr unsigned short m_nPhiSector = m_nPhi / m_phiGeo; 
 
    // Data structure/binning of the compressed hit pattern files (.txt.gz files)
    static constexpr unsigned short m_nPhiComp = 15; 
    static constexpr SectorBinning m_compAxialBins = {m_nOmega, m_nPhiComp, 1, m_nAxial, m_nPrio}; 
    static constexpr SectorBinning m_compStereoBins = {m_nOmega, m_nPhiComp, m_nCot, m_nStereo, m_nPrio}; 
 
    // Acceptance ranges + slot sizes to convert bins to track parameters (for getBinToVal method)
    static constexpr std::array<double, 2> m_omegaRange = {-4., 4.}; 
    static constexpr std::array<double, 2> m_phiRange = {0., 11.25}; 
    // These are optimized for minSuperAxial = 3 and minSuperStereo = 2
    static constexpr std::array<double, 2> m_cotRange = {2.3849627654510415, -1.0061730449796316}; 
    static constexpr double m_binSizeOmega = (m_omegaRange[1] - m_omegaRange[0]) / m_nOmega; 
    static constexpr double m_binSizePhi = (m_phiRange[1] - m_phiRange[0]) / m_nPhiSector; 
    static constexpr double m_binSizeCot = (m_cotRange[1] - m_cotRange[0]) / m_nCot; 
    static constexpr std::array<std::array<double, 2>, 3> m_acceptanceRanges = {m_omegaRange, m_phiRange, m_cotRange};
    static constexpr std::array<double, 3> m_binSizes = {m_binSizeOmega, m_binSizePhi, m_binSizeCot};
 
    // Array pointers to the hit representations and Hough space
    c2array* m_hitToSectorIDs = nullptr;
    c5array* m_compAxialHitReps = nullptr;
    c5array* m_compStereoHitReps = nullptr;
    c5array* m_expAxialHitReps = nullptr;
    c5array* m_expStereoHitReps = nullptr;
    c3array* m_houghSpace = nullptr;
  };

}