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 <cmath>
#include "trg/cdc/NDFinderDefs.h"
#include "trg/cdc/Clusterizend.h"
#include <Math/Vector3D.h>
 
namespace Belle2 {
  struct cellweight {
    cell_index index;
    ushort weight;
  };
 
  class NDFinderTrack {
  public:
    NDFinderTrack(std::vector<double> values, const SimpleCluster& cluster,
                  std::vector<ROOT::Math::XYZVector> houghspace, std::vector<ROOT::Math::XYZVector> ndreadout):
      m_cluster(cluster)
    {
      m_omega = values[0];
      m_phi = values[1];
      m_cotTheta = values[2];
      m_houghspace = houghspace;
      m_ndreadout = ndreadout;
    }
 
    virtual ~NDFinderTrack() {}
    double getOmega()
    {
      return m_omega;
    }
    double getPhi0()
    {
      return m_phi;
    }
    double getCot()
    {
      return m_cotTheta;
    }
    double getNRelHits()
    {
      return m_cluster.getHits().size();
    }
    std::vector<unsigned short> getRelHits()
    {
      return m_cluster.getHits();
    }
    std::vector<unsigned short> getRelHitsWeights()
    {
      return m_cluster.getWeights();
    }
    void reset()
    {
      m_omega = 0.;
      m_phi = 0.;
      m_cotTheta = 0.;
    }
    SimpleCluster getCluster()
    {
      return m_cluster;
    }
    std::vector<ROOT::Math::XYZVector> getHoughSpace()
    {
      return m_houghspace;
    }
    std::vector<ROOT::Math::XYZVector> getNDReadout()
    {
      return m_ndreadout;
    }
 
  private:
    double m_omega;
    double m_phi;
    double m_cotTheta;
    std::vector<unsigned short> m_relHits;
    std::vector<double> m_hitWeights;
    SimpleCluster m_cluster;
    std::vector<ROOT::Math::XYZVector> m_houghspace;
    std::vector<ROOT::Math::XYZVector> m_ndreadout;
  };
 
 
  class NDFinder {
  public:
 
    struct ndparameters {
      std::string axialFile = "data/trg/cdc/ndFinderArrayAxialComp.txt.gz";
      std::string stereoFile = "data/trg/cdc/ndFinderArrayStereoComp.txt.gz";
 
      unsigned char minSuperAxial  = 4;
 
      unsigned char minSuperStereo  = 3;
 
      float thresh = 0.85;
 
      unsigned char minCells  = 1;
 
      bool dbscanning = false;
 
      unsigned short minTotalWeight = 450;
 
      unsigned short minPeakWeight = 32;
 
      unsigned char iterations = 5;
 
      unsigned char omegaTrim = 5;
 
      unsigned char phiTrim = 4;
 
      unsigned char thetaTrim = 4;
 
      unsigned short phigeo = 32;
 
      unsigned short parcels = 7;
 
      unsigned short parcelsExp = 11;
    };
    std::vector<std::vector<float>> m_acceptRanges;
    std::vector<float> m_slotSizes;
 
    NDFinder() {}
 
    virtual ~NDFinder()
    {
      delete m_parrayAxial;
      delete m_parrayStereo;
      delete m_phoughPlane;
      delete m_parrayHitMod;
      delete m_pcompAxial;
      delete m_pcompStereo;
      delete m_parrayAxialExp;
      delete m_parrayStereoExp;
    }
 
 
    void init(unsigned short minWeight, unsigned char minPts, bool diagonal,
              unsigned char minSuperAxial, unsigned char minSuperStereo, float thresh,
              unsigned char minCells, bool dbscanning, unsigned short minTotalWeight, unsigned short minPeakWeight, unsigned char iterations,
              unsigned char omegaTrim, unsigned char phiTrim, unsigned char thetaTrim,
              bool verbose, std::string& axialFile, std::string& stereoFile);
 
    void initBins();
 
    void loadArray(const std::string& filename, ndbinning bins, c5array& hitsToTracks);
    void restoreZeros(ndbinning zerobins, ndbinning compbins, c5array& expArray, const c5array& compArray);
 
    void squeezeOne(c5array& writeArray, c5array& readArray, int outparcels, int inparcels, c5index ihit, c5index iprio, c5index itheta,
                    c5elem nomega);
 
    void squeezeAll(ndbinning writebins, c5array& writeArray, c5array& readArray, int outparcels, int inparcels);
 
    void initHitModAxSt(c2array& hitMod);
 
    void reset()
    {
      m_NDFinderTracks.clear();
      m_hitIds.clear();
      m_hitSLIds.clear();
      m_prioPos.clear();
      m_prioTime.clear();
      m_nHits = 0;
      m_vecDstart.clear();
      m_hitOrients.clear();
      delete m_phoughPlane;
      m_phoughPlane = new c3array(m_pc3shape);
    }
 
    void printParams();
 
    void addHit(unsigned short hitId, unsigned short hitSLId, unsigned short hitPrioPos, long hitPrioTime)
    {
      if (hitPrioPos > 0) { // skip "no hit"
        m_hitIds.push_back(hitId);
        m_hitSLIds.push_back(hitSLId);
        m_prioPos.push_back(3 - hitPrioPos);
        m_prioTime.push_back(hitPrioTime);
        m_nHits++;
      }
    }
 
    void findTracks();
 
    std::vector<NDFinderTrack>* getFinderTracks()
    {
      return &m_NDFinderTracks;
    }
 
    void printArray3D(c3array& hitsToTracks, ndbinning bins, ushort, ushort, ushort, ushort);
 
  protected:
 
    void getCM();
 
    void addLookup(unsigned short ihit);
 
    void addC3Comp(ushort hitr, ushort prio, const c5array& hitsToTracks,
                   short Dstart, ndbinning bins);
 
    std::vector<std::vector<unsigned short>> getHitsVsClusters(
                                            std::vector<SimpleCluster>& clusters);
 
    cell_index getMax(const std::vector<cell_index>&);
 
    ushort hitContrib(cell_index peak, ushort ihit);
 
    std::vector<SimpleCluster> allHitsToClusters(
      std::vector<std::vector<unsigned short>>& hitsVsClusters,
      std::vector<SimpleCluster>& clusters);
 
    std::vector<cellweight> getHighWeight(std::vector<cell_index> entries, float cutoff);
 
    std::vector<double> getWeightedMean(std::vector<cellweight>);
 
    std::vector<double> getBinToVal(std::vector<double>);
 
    double cdcTrackRadius(double pt)
    {
      return pt * 1e11 / (3e8 * 1.5); // div (c * B)
    }
 
    float transformVar(float estVal, int idx);
    std::vector<double> transform(std::vector<double> estimate);
 
  private:
 
    boost::array<c5index, 5> m_pc5shapeax;
    boost::array<c5index, 5> m_pc5shapest;
    boost::array<c3index, 3> m_pc3shape;
    boost::array<c2index, 2> m_pc2shapeHitMod;
    boost::array<c5index, 5> m_pc5shapeCompAx;
    boost::array<c5index, 5> m_pc5shapeCompSt;
    boost::array<c5index, 5> m_pc5shapeExpAx;
    boost::array<c5index, 5> m_pc5shapeExpSt;
 
    c5array* m_parrayAxial = nullptr;
    c5array* m_parrayStereo = nullptr;
    c3array* m_phoughPlane = nullptr;
    c2array* m_parrayHitMod = nullptr;
    c5array* m_pcompAxial = nullptr;
    c5array* m_pcompStereo = nullptr;
    c5array* m_parrayAxialExp = nullptr;
    c5array* m_parrayStereoExp = nullptr;
 
    std::vector<int> m_nWires;
 
    std::vector<NDFinderTrack> m_NDFinderTracks;
 
    std::vector<unsigned short> m_hitIds;
 
    std::vector<unsigned short> m_hitSLIds;
 
    std::vector<unsigned short> m_prioPos;
 
    std::vector<long> m_prioTime;
 
    std::vector<unsigned short> m_hitOrients;
 
    std::vector<short> m_vecDstart;
 
    unsigned short m_nHits{0};
 
    ndparameters m_params;
 
    ndbinning m_axBins;
    ndbinning m_stBins;
    ndbinning m_fullBins;
    ndbinning m_compAxBins;
    ndbinning m_compStBins;
    ndbinning m_expAxBins;
    ndbinning m_expStBins;
 
    clustererParams m_clustererParams;
    std::vector<ushort> m_planeShape;
 
    unsigned short m_nPhiFull{0};
    unsigned short m_nPhiOne{0};
    unsigned short m_nPhiComp{0};
    unsigned short m_nPhiExp{0};
    unsigned short m_nPhiUse{0};
    unsigned short m_nOmega{0};
    unsigned short m_nTheta{0};
 
    unsigned short m_nSL{0};
    unsigned short m_nTS{0};
    unsigned short m_nAx{0};
    unsigned short m_nSt{0};
    unsigned short m_nPrio{0};
 
    Belle2::Clusterizend m_clusterer;
    Belle2::Clusterizend m_clusterer2;
 
    bool m_verbose{false};
  };
}