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.                  *
 **************************************************************************/
 
#ifndef NDFINDER_H
#define NDFINDER_H
 
#include <cmath>
#include "trg/cdc/NDFinderDefs.h"
#include "trg/cdc/Clusterizend.h"
 
namespace Belle2 {
  struct cellweight {
    cell_index index;
    ushort weight;
  };
 
  class NDFinderTrack {
  public:
    NDFinderTrack(std::vector<double> values, const SimpleCluster& cluster):
      m_cluster(cluster)
    {
      m_omega = values[0];
      m_phi = values[1];
      m_cotTheta = values[2];
    }
 
    virtual ~NDFinderTrack() {}
    double getOmega()
    {
      return m_omega;
    }
    double getPhi0()
    {
      return m_phi;
    }
    double getCot()
    {
      return m_cotTheta;
    }
    double get_nRelHits()
    {
      return m_cluster.get_hits().size();
    }
    std::vector<unsigned short> get_relHits()
    {
      return m_cluster.get_hits();
    }
    std::vector<unsigned short> get_relHitsWeights()
    {
      return m_cluster.get_weights();
    }
    void reset()
    {
      m_omega = 0.;
      m_phi = 0.;
      m_cotTheta = 0.;
    }
 
  private:
    double m_omega;
    double m_phi;
    double m_cotTheta;
    std::vector<unsigned short> m_relHits;
    std::vector<double> m_hitWeights;
    SimpleCluster m_cluster;
  };
 
 
  class NDFinder {
  public:
 
    struct ndparameters {
      std::string axialFile = "data/trg/cdc/ndFinderArrayAxialComp.txt.gz";
      std::string stereoFile = "data/trg/cdc/ndFinderArrayStereoComp.txt.gz";
 
      long unsigned int minhits  = 4;
      long unsigned int minhits_axial  = 0;
      float thresh            = 0.85;
 
      float minassign = 0.2;
 
      long unsigned int mincells  = 5;
 
      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(int minweight, int minpts, bool diagonal,
              int minhits, int minhits_axial, double thresh, double minassign,
              int mincells, 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_prioPos.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 hitPrioPos)
    {
      if (hitPrioPos > 0) { // skip "no hit"
        m_hitIds.push_back(hitId);
        m_prioPos.push_back(3 - hitPrioPos);
        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> relateHitsToClusters(
      std::vector<std::vector<unsigned short>>& hitsVsClusters,
      std::vector<SimpleCluster>& clusters);
 
    int hitToCluster(std::vector<std::vector<unsigned short>>& hitsVsClusters,
                     unsigned short hit);
 
    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 cdc_track_radius(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_prioPos;
 
    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;
 
    clusterer_params m_clusterer_params;
    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;
 
    bool m_verbose{false};
  };
}
#endif