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/ActsGeometry.h>
0014 #include <trackbase/TrkrClusterCrossingAssoc.h>
0015 #include <trackbase/TrkrDefs.h>  // for cluskey
0016
0017 #include <g4detectors/PHG4CylinderGeom.h>
0018 #include <g4detectors/PHG4CylinderGeomContainer.h>
0019
0020 #include <g4mvtx/PHG4MvtxMisalignment.h>
0021
0022 #pragma GCC diagnostic push
0023 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0024 #include <boost/geometry/geometries/box.hpp>  // for box
0025 #pragma GCC diagnostic pop
0026
0027 #include <boost/geometry/geometries/point.hpp>  // for point
0028 #include <boost/geometry/index/rtree.hpp>       // for ca
0029
0030 #include <cmath>    // for M_PI
0031 #include <cstdint>  // for uint64_t
0032 #include <map>      // for map
0033 #include <memory>
0034 #include <set>
0035 #include <string>  // for string
0036 #include <unordered_map>
0037 #include <utility>  // for pair
0038 #include <vector>   // for vector
0039
0040 class PHCompositeNode;
0041 class TrkrCluster;
0042 class TrackSeed;
0043 class TrackSeed_v2;
0044
0045 namespace bg = boost::geometry;
0046 namespace bgi = boost::geometry::index;
0047
0048 class PHCASiliconSeeding : public PHTrackSeeding
0049 {
0050  public:
0051   using point = bg::model::point<float, 2, bg::cs::cartesian>;
0052   using box = bg::model::box<point>;
0053   using pointKey = std::pair<point, TrkrDefs::cluskey>;                 // phi and z in the key
0054   using coordKey = std::pair<std::array<float, 2>, TrkrDefs::cluskey>;  // just use phi and Z, no longer needs the layer
0055
0056   using keyList = std::vector<TrkrDefs::cluskey>;
0057   using keyLists = std::vector<keyList>;
0058   using keyListPerLayer = std::array<keyList, 55>;
0059
0060   using PositionMap = std::unordered_map<TrkrDefs::cluskey, Acts::Vector3>;
0061
0062   PHCASiliconSeeding(
0063       const std::string& name = "PHCASiliconSeeding",
0064       unsigned int start_layer = 7,
0065       unsigned int end_layer = 55,
0066       unsigned int min_clusters_per_track = 5,
0067       float neighbor_phi_width = .02,
0068       float eta_allowance = 1.1);
0069
0070   ~PHCASiliconSeeding() override {}
0071   void SetLayerRange(unsigned int layer_low, unsigned int layer_up)
0072   {
0073     _start_layer = layer_low;
0074     _end_layer = layer_up;
0075   }
0076   void SetMVTXStrobeIDRange(int low_strobe, int high_strobe)
0077   {
0078     _lowest_allowed_strobeid = low_strobe;
0079     _highest_allowed_strobeid = high_strobe;
0080   }
0081   void SetSearchWindow(float eta_allowance, float phi_width)
0082   {
0083     _eta_allowance = eta_allowance;
0084     _drdz_allowance = (eta_allowance < 1E-5) ? (2. * std::exp(-5.0)) / (1 - std::exp(-1. * 5.0 * 2)) : (2. * std::exp(-eta_allowance)) / (1 - std::exp(-1. * eta_allowance * 2));  // If eta_allowance is very small, set to the value for eta = 5.0 for large allowance
0085     _neighbor_phi_width = phi_width;
0086   }
0087   void SetPropagateMaxDCAxy(float dcaxy)
0088   {
0089     _propagate_max_dcaxy = dcaxy;
0090   }
0091   void SetPropagateMaxDCAz(float dcaz)
0092   {
0093     _propagate_max_dcaz = dcaz;
0094   }
0095   void SetMaxCosTripletBreakingAngle(float cos_angle)
0096   {
0097     _max_cos_angle = cos_angle;
0098   }
0099   void SetAlgoUseBestTriplet(bool use_best)
0100   {
0101     _use_best = use_best;
0102   }
0103   void SetRequireINTTConsistency(bool req)
0104   {
0105     _require_INTT_consistency = req;
0106   }
0107   void SetMinClustersPerTrack(unsigned int minClus)
0108   {
0109     _min_clusters_per_track = minClus;
0110   }
0111   void SetMinMVTXClusters(unsigned int minMVTX)
0112   {
0113     _min_mvtx_clusters = minMVTX;
0114   }
0115   void SetMinINTTClusters(unsigned int minINTT)
0116   {
0117     _min_intt_clusters = minINTT;
0118   }
0119   void SetTrackMapName(const std::string& trackmap_name)
0120   {
0121     _module_trackmap_name = trackmap_name;
0122   }
0123   void Identify() const
0124   {
0125     std::cout << "----- Configuration parameters -----" << std::endl;
0126     std::cout << " - Start layer: " << _start_layer << std::endl;
0127     std::cout << " - End layer: " << _end_layer << std::endl;
0128     std::cout << " - MVTX strobe ID range: [" << _lowest_allowed_strobeid << ", " << _highest_allowed_strobeid << "]" << std::endl;
0129     std::cout << " - Search window eta allowance: " << _eta_allowance << ", drdz allowance: " << _drdz_allowance << std::endl;
0130     std::cout << " - Search window neighbor phi width: " << _neighbor_phi_width << std::endl;
0131     std::cout << " - Propagate max DCAxy: " << _propagate_max_dcaxy << std::endl;
0132     std::cout << " - Propagate max DCAz: " << _propagate_max_dcaz << std::endl;
0133     std::cout << " - Max cos(triplet breaking angle): " << _max_cos_angle << std::endl;
0134     std::cout << " - Use best triplet: " << (_use_best ? "true" : "false") << std::endl;
0135     std::cout << " - Require INTT consistency: " << (_require_INTT_consistency ? "true" : "false") << std::endl;
0136     std::cout << " - Minimum clusters per track: " << _min_clusters_per_track << std::endl;
0137     std::cout << " - Minimum MVTX clusters: " << _min_mvtx_clusters << std::endl;
0138     std::cout << " - Minimum INTT clusters: " << _min_intt_clusters << std::endl;
0139     std::cout << " - Track Map Name: " << _module_trackmap_name << std::endl;
0140     std::cout << "------------------------------------" << std::endl;
0141   }
0142
0143  protected:
0144   int Setup(PHCompositeNode* topNode) override;
0145   int Process(PHCompositeNode* topNode) override;
0146   int InitializeGeometry(PHCompositeNode* topNode);
0147   int End() override;
0148
0149  private:
0150   unsigned int _start_layer;
0151   unsigned int _end_layer;
0152   float _strobe_width = 9.9; // mus
0153   unsigned int _min_clusters_per_track;
0154   unsigned int _min_mvtx_clusters = 2;
0155   unsigned int _min_intt_clusters = 1;
0156
0157   float _eta_allowance = 1.1;
0158   float _drdz_allowance = (2. * std::exp(-_eta_allowance)) / (1 - std::exp(-1. * _eta_allowance * 2));  // default for eta allowance of 1.1
0159   float _neighbor_phi_width;
0160
0161   int _lowest_allowed_strobeid = -5;
0162   int _highest_allowed_strobeid = 5;
0163
0164   float _propagate_max_dcaxy = 0.4;
0165   float _propagate_max_dcaz = 1.;
0166
0167   float _max_cos_angle = -0.95;
0168   bool _use_best = true;
0169   bool _require_INTT_consistency = true;
0170
0171   std::array<float, 55> dphi_per_layer{};
0172   std::array<float, 55> max_dcaxy_perlayer{};
0173   std::array<float, 55> max_dcaz_perlayer{};
0174   std::array<float, 7> radius_per_layer{};  // radius of each layer
0175
0176   struct Triplet
0177   {
0178     TrkrDefs::cluskey bottom;
0179     TrkrDefs::cluskey center;
0180     TrkrDefs::cluskey top;
0181   };
0182
0183   /// get global position for a given cluster
0184   /**
0185    * uses ActsTransformation to convert cluster local position into global coordinates
0186    * incorporates TPC distortion correction, if present
0187    */
0188   bool timingMismatch(const TrackSeed& seed) const;
0189   Acts::Vector3 getGlobalPosition(TrkrDefs::cluskey, TrkrCluster*) const;
0190   std::pair<PositionMap, keyListPerLayer> FillGlobalPositions();
0191
0192   std::vector<std::vector<Triplet>> CreateLinks(const PHCASiliconSeeding::PositionMap& globalPositions, const PHCASiliconSeeding::keyListPerLayer& ckeys);
0193   std::vector<keyList> FollowLinks(const std::vector<std::vector<Triplet>>& triplets);
0194   std::vector<coordKey> FillTree(bgi::rtree<pointKey, bgi::quadratic<16>>&, const keyList&, const PositionMap&, int layer);
0195   int FindSeeds(const PositionMap&, const keyListPerLayer&);
0196
0197   void QueryTree(const bgi::rtree<pointKey, bgi::quadratic<16>>& rtree, double phimin, double zmin, double phimax, double zmax, std::vector<pointKey>& returned_values) const;
0198   float getSeedQuality(const TrackSeed_v2& seed, const PositionMap& globalPositions) const;
0199   void HelixPropagate(std::vector<TrackSeed_v2>& seeds, const PositionMap& globalPositions) const;
0200   void FitSeed(TrackSeed_v2& seed, const PositionMap& globalPositions) const;
0201   std::vector<TrackSeed_v2> FitSeeds(const std::vector<keyList>& seeds, const PositionMap& globalPositions) const;
0202
0203   std::set<short> GetINTTClusterCrossings(const TrkrDefs::cluskey ckey) const;
0204   short GetCleanINTTClusterCrossing(const TrkrDefs::cluskey ckey) const;
0205   int GetClusterTimeIndex(const TrkrDefs::cluskey ckey) const;
0206   bool ClusterTimesAreCompatible(const uint8_t trkr_id, const int time_index, const TrkrDefs::cluskey ckey) const;
0207   bool ClusterTimesAreCompatible(const TrkrDefs::cluskey clus_a, const TrkrDefs::cluskey clus_b) const;
0208
0209   void publishSeeds(const std::vector<TrackSeed_v2>& seeds) const;
0210
0211   // set up layer radii
0212   void SetupDefaultLayerRadius();
0213   float GetMvtxRadiusByPhi(float clusphi, int layer) const;
0214
0215   /// acts geometry
0216   ActsGeometry* m_tGeometry{nullptr};
0217   // cylinder geometry for mvtx and intt
0218   PHG4CylinderGeomContainer* geom_container_mvtx = nullptr;
0219   PHG4CylinderGeomContainer* geom_container_intt = nullptr;
0220   PHG4MvtxMisalignment* mvtxmisalignment = nullptr;
0221   std::vector<double> v_globaldisplacement = {0., 0., 0.};
0222   float radius_displacement = 0.;
0223
0224   //  TrackSeedContainer *m_seedContainer = nullptr;
0225   TrkrClusterContainer* m_clusterMap = nullptr;
0226   TrkrClusterCrossingAssoc* m_clusterCrossingMap = nullptr;
0227
0228   std::string _module_trackmap_name = "SiliconTrackSeedContainer";
0229
0230   std::vector<bgi::rtree<pointKey, bgi::quadratic<16>>> _rtrees;  // need three layers at a time
0231 };
0232
0233 #endif