packages/trackreco/PHCASiliconSeeding.h

0001 #ifndef TRACKRECO_PHCASILICONSEEDING_H
0002 #define TRACKRECO_PHCASILICONSEEDING_H
0003
0004 /*!
0005  *  \file PHCASiliconSeeding.cc
0006  *  \brief Silicon track seeding using ALICE-style "cellular automaton" (CA) algorithm
0007  *  \detail
0008  *  \author Michael Peters
0009  */
0010
0011 #include "PHTrackSeeding.h"  // for PHTrackSeeding
0012
0013 #include <trackbase/TrkrDefs.h>  // for cluskey
0014 #include <trackbase/TrkrClusterCrossingAssoc.h>
0015 #include <trackbase/ActsGeometry.h>
0016
0017 #pragma GCC diagnostic push
0018 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0019 #include <boost/geometry/geometries/box.hpp>  // for box
0020 #pragma GCC diagnostic pop
0021
0022 #include <boost/geometry/geometries/point.hpp>  // for point
0023 #include <boost/geometry/index/rtree.hpp>       // for ca
0024
0025 #include <cmath>    // for M_PI
0026 #include <cstdint>  // for uint64_t
0027 #include <map>      // for map
0028 #include <unordered_map>
0029 #include <memory>
0030 #include <set>
0031 #include <string>         // for string
0032 #include <utility>  // for pair
0033 #include <vector>   // for vector
0034
0035 class PHCompositeNode;
0036 class TrkrCluster;
0037 class TrackSeed_v2;
0038
0039 namespace bg = boost::geometry;
0040 namespace bgi = boost::geometry::index;
0041
0042 class PHCASiliconSeeding : public PHTrackSeeding
0043 {
0044  public:
0045   using point = bg::model::point<float, 2, bg::cs::cartesian>;
0046   using box = bg::model::box<point>;
0047   using pointKey = std::pair<point, TrkrDefs::cluskey>;                 // phi and z in the key
0048   using coordKey = std::pair<std::array<float, 2>, TrkrDefs::cluskey>;  // just use phi and Z, no longer needs the layer
0049
0050   using keyList = std::vector<TrkrDefs::cluskey>;
0051   using keyLists = std::vector<keyList>;
0052   using keyListPerLayer = std::array<keyList, 55>;
0053
0054   using PositionMap = std::unordered_map<TrkrDefs::cluskey, Acts::Vector3>;
0055
0056   PHCASiliconSeeding(
0057       const std::string& name = "PHCASiliconSeeding",
0058       unsigned int start_layer = 7,
0059       unsigned int end_layer = 55,
0060       unsigned int min_clusters_per_track = 5,
0061       float neighbor_phi_width = .02,
0062       float neighbor_z_width = .01
0063   );
0064
0065   ~PHCASiliconSeeding() override {}
0066   void SetLayerRange(unsigned int layer_low, unsigned int layer_up)
0067   {
0068     _start_layer = layer_low;
0069     _end_layer = layer_up;
0070   }
0071   void SetMVTXStrobeIDRange(int low_strobe, int high_strobe)
0072   {
0073     _lowest_allowed_strobeid = low_strobe;
0074     _highest_allowed_strobeid = high_strobe;
0075   }
0076   void SetSearchWindow(float z_width, float phi_width)
0077   {
0078     _neighbor_z_width = z_width;
0079     _neighbor_phi_width = phi_width;
0080   }
0081   void SetPropagateMaxDCAxy(float dcaxy)
0082   {
0083     _propagate_max_dcaxy = dcaxy;
0084   }
0085   void SetPropagateMaxDCAz(float dcaz)
0086   {
0087     _propagate_max_dcaz = dcaz;
0088   }
0089   void SetMaxCosTripletBreakingAngle(float cos_angle)
0090   {
0091     _max_cos_angle = cos_angle;
0092   }
0093   void SetAlgoUseBestTriplet(bool use_best)
0094   {
0095     _use_best = use_best;
0096   }
0097   void SetRequireINTTConsistency(bool req)
0098   {
0099     _require_INTT_consistency = req;
0100   }
0101   void SetMinClustersPerTrack(unsigned int minClus)
0102   {
0103     _min_clusters_per_track = minClus;
0104   }
0105
0106  protected:
0107   int Setup(PHCompositeNode* topNode) override;
0108   int Process(PHCompositeNode* topNode) override;
0109   int InitializeGeometry(PHCompositeNode* topNode);
0110   int End() override;
0111
0112  private:
0113   unsigned int _start_layer;
0114   unsigned int _end_layer;
0115   unsigned int _min_clusters_per_track;
0116   float _neighbor_phi_width;
0117   float _neighbor_z_width;
0118
0119   int _lowest_allowed_strobeid = -5;
0120   int _highest_allowed_strobeid = 5;
0121
0122   float _propagate_max_dcaxy = 0.4;
0123   float _propagate_max_dcaz = 1.;
0124
0125   float _max_cos_angle = -0.95;
0126   bool _use_best = true;
0127   bool _require_INTT_consistency = true;
0128
0129   std::array<float, 55> dZ_per_layer{};
0130   std::array<float, 55> dphi_per_layer{};
0131   std::array<float, 55> max_dcaxy_perlayer{};
0132   std::array<float, 55> max_dcaz_perlayer{};
0133
0134   struct Triplet
0135   {
0136     TrkrDefs::cluskey bottom;
0137     TrkrDefs::cluskey center;
0138     TrkrDefs::cluskey top;
0139   };
0140
0141   /// get global position for a given cluster
0142   /**
0143    * uses ActsTransformation to convert cluster local position into global coordinates
0144    * incorporates TPC distortion correction, if present
0145    */
0146   Acts::Vector3 getGlobalPosition(TrkrDefs::cluskey, TrkrCluster*) const;
0147   std::pair<PositionMap, keyListPerLayer> FillGlobalPositions();
0148
0149   std::vector<std::vector<Triplet>> CreateLinks(const PHCASiliconSeeding::PositionMap& globalPositions, const PHCASiliconSeeding::keyListPerLayer& ckeys);
0150   std::vector<keyList> FollowLinks(const std::vector<std::vector<Triplet>>& triplets);
0151   std::vector<coordKey> FillTree(bgi::rtree<pointKey, bgi::quadratic<16>>&, const keyList&, const PositionMap&, int layer);
0152   int FindSeeds(const PositionMap&, const keyListPerLayer&);
0153
0154   void QueryTree(const bgi::rtree<pointKey, bgi::quadratic<16>>& rtree, double phimin, double zmin, double phimax, double zmax, std::vector<pointKey>& returned_values) const;
0155   float getSeedQuality(const TrackSeed_v2& seed, const PositionMap& globalPositions) const;
0156   void HelixPropagate(std::vector<TrackSeed_v2>& seeds, const PositionMap& globalPositions) const;
0157   void FitSeed(TrackSeed_v2& seed, const PositionMap& globalPositions) const;
0158   std::vector<TrackSeed_v2> FitSeeds(const std::vector<keyList>& seeds, const PositionMap& globalPositions) const;
0159
0160   std::set<short> GetINTTClusterCrossings(const TrkrDefs::cluskey ckey) const;
0161   short GetCleanINTTClusterCrossing(const TrkrDefs::cluskey ckey) const;
0162   int GetClusterTimeIndex(const TrkrDefs::cluskey ckey) const;
0163   bool ClusterTimesAreCompatible(const uint8_t trkr_id, const int time_index, const TrkrDefs::cluskey ckey) const;
0164   bool ClusterTimesAreCompatible(const TrkrDefs::cluskey clus_a, const TrkrDefs::cluskey clus_b) const;
0165
0166   void publishSeeds(const std::vector<TrackSeed_v2>& seeds) const;
0167
0168   /// acts geometry
0169   ActsGeometry* m_tGeometry{nullptr};
0170
0171   std::string trackmapname = "SiliconTrackSeedContainer";
0172   //  TrackSeedContainer *m_seedContainer = nullptr;
0173   TrkrClusterContainer *m_clusterMap = nullptr;
0174   TrkrClusterCrossingAssoc *m_clusterCrossingMap = nullptr;
0175
0176   std::vector<bgi::rtree<pointKey, bgi::quadratic<16>>> _rtrees;  // need three layers at a time
0177 };
0178
0179 #endif