packages/trackreco/PHCASeeding.h

0001 #ifndef TRACKRECO_PHCASEEDING_H
0002 #define TRACKRECO_PHCASEEDING_H
0003
0004 /*!
0005  *  \file PHCASeeding.cc
0006  *  \brief Track seeding using ALICE-style "cellular automaton" (CA) algorithm
0007  *  \detail
0008  *  \author Michael Peters & Christof Roland
0009  */
0010
0011 // Statements for if we want to save out the intermediary clustering steps
0012 /* #define _PHCASEEDING_CLUSTERLOG_TUPOUT_ */
0013 /* #define _PHCASEEDING_CHAIN_FORKS_ */
0014 /* #define _PHCASEEDING_TIMER_OUT_ */
0015
0016 #include "PHTrackSeeding.h"  // for PHTrackSeeding
0017
0018 #include <tpc/TpcGlobalPositionWrapper.h>
0019
0020 #include <trackbase/TrkrDefs.h>  // for cluskey
0021 #include <trackbase_historic/TrackSeed_v2.h>
0022
0023 #include <phool/PHTimer.h>  // for PHTimer
0024
0025 #include <Eigen/Core>
0026 #include <Eigen/Dense>
0027
0028 #pragma GCC diagnostic push
0029 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0030 #include <boost/geometry/geometries/box.hpp>  // for box
0031 #pragma GCC diagnostic pop
0032
0033 #include <boost/geometry/geometries/point.hpp>  // for point
0034 #include <boost/geometry/index/rtree.hpp>       // for ca
0035
0036 #include <cmath>    // for M_PI
0037 #include <cstdint>  // for uint64_t
0038 #include <map>      // for map
0039 #include <memory>
0040 #include <set>
0041 #include <string>         // for string
0042 #include <unordered_map>  // for map
0043 #include <unordered_set>
0044 #include <utility>  // for pair
0045 #include <vector>   // for vector
0046
0047 #include <TFile.h>
0048 #include <TNtuple.h>
0049
0050 class ActsGeometry;
0051 class PHCompositeNode;
0052 class PHTimer;
0053 class SvtxTrack_v3;
0054 class TpcDistortionCorrectionContainer;
0055 class TrkrCluster;
0056
0057 namespace bg = boost::geometry;
0058 namespace bgi = boost::geometry::index;
0059
0060 class PHCASeeding : public PHTrackSeeding
0061 {
0062  public:
0063   static const int _NLAYERS_TPC = 48;
0064   static const int _FIRST_LAYER_TPC = 7;
0065   // move `using` statements inside of the class to avoid polluting the global namespace
0066
0067   using point = bg::model::point<float, 2, bg::cs::cartesian>;
0068   using box = bg::model::box<point>;
0069   using pointKey = std::pair<point, TrkrDefs::cluskey>;                 // phi and z in the key
0070   using coordKey = std::pair<std::array<float, 2>, TrkrDefs::cluskey>;  // just use phi and Z, no longer needs the layer
0071
0072   using keyList = std::vector<TrkrDefs::cluskey>;
0073   using keyLists = std::vector<keyList>;
0074   using keyListPerLayer = std::array<keyList, _NLAYERS_TPC>;
0075   using keySet = std::set<TrkrDefs::cluskey>;
0076
0077   using keyLink = std::pair<TrkrDefs::cluskey, TrkrDefs::cluskey>;
0078   using keyLinks = std::vector<keyLink>;
0079   using keyLinkPerLayer = std::array<std::vector<keyLink>, _NLAYERS_TPC>;
0080
0081   using PositionMap = std::unordered_map<TrkrDefs::cluskey, Acts::Vector3>;
0082
0083   std::array<float, 55> dZ_per_layer{};
0084   std::array<float, 55> dphi_per_layer{};
0085
0086   PHCASeeding(
0087       const std::string& name = "PHCASeeding",
0088       unsigned int start_layer = 7,
0089       unsigned int end_layer = 55,
0090       unsigned int min_nhits_per_cluster = 0,
0091       unsigned int min_clusters_per_track = 5,
0092       const unsigned int nlayers_maps = 3,
0093       const unsigned int nlayers_intt = 4,
0094       const unsigned int nlayers_tpc = 48,
0095       float neighbor_phi_width = .02,
0096       float neighbor_z_width = .01,
0097       float maxSinPhi = 0.999
0098       /* float cosTheta_limit = -0.8 */
0099   );
0100
0101   ~PHCASeeding() override {}
0102   void SetSplitSeeds(bool opt = true) { _split_seeds = opt; }
0103   void SetLayerRange(unsigned int layer_low, unsigned int layer_up)
0104   {
0105     _start_layer = layer_low;
0106     _end_layer = layer_up;
0107   }
0108   void SetSearchWindow(float z_width, float phi_width)
0109   {
0110     _neighbor_z_width = z_width;
0111     _neighbor_phi_width = phi_width;
0112   }
0113   void SetClusAdd_delta_window(float _dzdr_cutoff, float _dphidr2_cutoff)
0114   {
0115     _clusadd_delta_dzdr_window = _dzdr_cutoff;
0116     _clusadd_delta_dphidr2_window = _dphidr2_cutoff;
0117   }
0118   void SetMinHitsPerCluster(unsigned int minHits) { _min_nhits_per_cluster = minHits; }
0119   void SetMinClustersPerTrack(unsigned int minClus) { _min_clusters_per_track = minClus; }
0120   void SetNClustersPerSeedRange(unsigned int minClus, unsigned int maxClus)
0121   {
0122     _min_clusters_per_seed = minClus;
0123     _max_clusters_per_seed = maxClus;
0124     if (_min_clusters_per_seed < 3)
0125     {
0126       std::cout << " Error in SetNClustersPerSeedRange: " << __FILE__
0127                 << " min value cannot be less than three." << std::endl;
0128       assert(false);
0129     }
0130   }
0131
0132   // obsolete
0133   void set_field_dir(const double) {}
0134   void magFieldFile(const std::string&) {}
0135   void useConstBField(bool) {}
0136   void constBField(float) {}
0137
0138   void useFixedClusterError(bool opt) { _use_fixed_clus_err = opt; }
0139   void setFixedClusterError(int i, double val) { _fixed_clus_err.at(i) = val; }
0140   void set_pp_mode(bool mode) { _pp_mode = mode; }
0141   void reject_zsize1_clusters(bool mode){_reject_zsize1 = mode;}
0142   void setNeonFraction(double frac) { Ne_frac = frac; };
0143   void setArgonFraction(double frac) { Ar_frac = frac; };
0144   void setCF4Fraction(double frac) { CF4_frac = frac; };
0145   void setNitrogenFraction(double frac) { N2_frac = frac; };
0146   void setIsobutaneFraction(double frac) { isobutane_frac = frac; };
0147
0148  protected:
0149   int Setup(PHCompositeNode* topNode) override;
0150   int Process(PHCompositeNode* topNode) override;
0151   int InitializeGeometry(PHCompositeNode* topNode);
0152   /* int FindSeedsLayerSkip(double cosTheta_limit); */
0153   int End() override;
0154
0155  private:
0156   bool _save_clus_proc = false;
0157   TFile* _f_clustering_process = nullptr;
0158   int _tupout_count = -1;
0159   int _n_tupchains = -1;
0160   // Save the steps of the clustering
0161   TNtuple* _tupclus_all = nullptr;          // all clusters
0162   TNtuple* _tupclus_links = nullptr;        //  clusters which are linked
0163   TNtuple* _tupclus_bilinks = nullptr;      // bi-linked clusters
0164   TNtuple* _tupclus_seeds = nullptr;        // seed (outermost three bi-linked chain
0165   TNtuple* _tupclus_grown_seeds = nullptr;  // seeds saved out
0166                                             //
0167   TNtuple* _tup_chainbody = nullptr;
0168   TNtuple* _tup_chainfork = nullptr;
0169
0170   TNtuple* _tupwin_link = nullptr;       // cluster pairs with dphi and dz (very wide windows) x0,x1,y0,y1,z0,z1
0171   TNtuple* _tupwin_cos_angle = nullptr;  // directed cosine (all cluster triples which are found) x0,x1,x2
0172   TNtuple* _tupwin_seed23 = nullptr;     // from third to fourth link -- can only use passing seeds
0173   TNtuple* _tupwin_seedL1 = nullptr;     // from third to fourth link -- can only use passing seeds
0174
0175   TNtuple* _search_windows = nullptr;  // This is really just a lazy way to store what search paramaters where used.
0176                                        // It would be equally valid in a TMap or map<string,float>
0177   // functions used to fill tuples -- only defined if _PHCASEEDING_CLUSTERLOG_TUPOUT_ is defined in preprocessor
0178   void write_tuples();
0179   void fill_tuple(TNtuple*, float, TrkrDefs::cluskey, const Acts::Vector3&) const;
0180   void fill_tuple_with_seed(TNtuple*, const keyList&, const PositionMap&) const;
0181   void process_tupout_count();
0182   void FillTupWinLink(bgi::rtree<pointKey, bgi::quadratic<16>>&, const coordKey&, const PositionMap&) const;
0183   void FillTupWinCosAngle(const TrkrDefs::cluskey, const TrkrDefs::cluskey, const TrkrDefs::cluskey, const PositionMap&, double cos_angle, bool isneg) const;
0184   void FillTupWinGrowSeed(const keyList& seed, const keyLink& link, const PositionMap& globalPositions) const;
0185   void fill_split_chains(const keyList& chain, const keyList& keylinks, const PositionMap& globalPositions, int& nchains) const;
0186   /* void fill_tuple_with_seed(TN */
0187
0188   // have a comparator to search vector of sorted bilinks for links with starting
0189   // a key of a given value
0190   class CompKeyToBilink
0191   {
0192    public:
0193     bool operator()(const keyLink& p, const TrkrDefs::cluskey& val) const
0194     {
0195       return p.first < val;
0196     }
0197     bool operator()(const TrkrDefs::cluskey& val, const keyLink& p) const
0198     {
0199       return val < p.first;
0200     }
0201   };
0202   std::pair<std::vector<keyLink>::iterator, std::vector<keyLink>::iterator> FindBilinks(const TrkrDefs::cluskey& key);
0203
0204   /// tpc distortion correction utility class
0205   TpcDistortionCorrection m_distortionCorrection;
0206
0207   /// get global position for a given cluster
0208   /**
0209    * uses ActsTransformation to convert cluster local position into global coordinates
0210    * incorporates TPC distortion correction, if present
0211    */
0212   Acts::Vector3 getGlobalPosition(TrkrDefs::cluskey, TrkrCluster*) const;
0213   std::pair<PositionMap, keyListPerLayer> FillGlobalPositions();
0214   std::pair<keyLinks, keyLinkPerLayer> CreateBiLinks(const PositionMap& globalPositions, const keyListPerLayer& ckeys);
0215   PHCASeeding::keyLists FollowBiLinks(const keyLinks& trackSeedPairs, const keyLinkPerLayer& bilinks, const PositionMap& globalPositions) const;
0216   std::vector<coordKey> FillTree(bgi::rtree<pointKey, bgi::quadratic<16>>&, const keyList&, const PositionMap&, int layer);
0217   int FindSeedsWithMerger(const PositionMap&, const keyListPerLayer&);
0218
0219   void QueryTree(const bgi::rtree<pointKey, bgi::quadratic<16>>& rtree, double phimin, double zmin, double phimax, double zmax, std::vector<pointKey>& returned_values) const;
0220   std::vector<TrackSeed_v2> RemoveBadClusters(const std::vector<keyList>& seeds, const PositionMap& globalPositions) const;
0221   double getMengerCurvature(TrkrDefs::cluskey a, TrkrDefs::cluskey b, TrkrDefs::cluskey c, const PositionMap& globalPositions) const;
0222
0223   void publishSeeds(const std::vector<TrackSeed_v2>& seeds) const;
0224
0225   // int _nlayers_all;
0226   // unsigned int _nlayers_seeding;
0227   // std::vector<int> _seeding_layer;
0228
0229   const unsigned int _nlayers_maps;
0230   const unsigned int _nlayers_intt;
0231   const unsigned int _nlayers_tpc;
0232   unsigned int _start_layer;
0233   unsigned int _end_layer;
0234   unsigned int _min_nhits_per_cluster;
0235   unsigned int _min_clusters_per_track;
0236   unsigned int _max_clusters_per_seed = 6;  // currently not used
0237   unsigned int _min_clusters_per_seed = 6;  // currently the abs. number used
0238   //  float _cluster_z_error;
0239   //  float _cluster_alice_y_error;
0240   float _neighbor_phi_width;
0241   float _neighbor_z_width;
0242   float _clusadd_delta_dzdr_window = 0.5;
0243   float _clusadd_delta_dphidr2_window = 0.005;
0244   float _max_sin_phi;
0245   /* float _cosTheta_limit; */
0246   double _rz_outlier_threshold = 0.1;
0247   double _xy_outlier_threshold = 0.1;
0248   bool _split_seeds = true;
0249   bool _reject_zsize1 = false;
0250   bool _use_fixed_clus_err = false;
0251   bool _pp_mode = false;
0252   std::array<double, 3> _fixed_clus_err = {.1, .1, .1};
0253
0254   /// acts geometry
0255   ActsGeometry* m_tGeometry{nullptr};
0256
0257   /// global position wrapper
0258   TpcGlobalPositionWrapper m_globalPositionWrapper;
0259
0260   std::unique_ptr<PHTimer> t_seed;
0261   std::unique_ptr<PHTimer> t_fill;
0262   std::unique_ptr<PHTimer> t_makebilinks;
0263   std::unique_ptr<PHTimer> t_makeseeds;
0264   /* std::array<bgi::rtree<pointKey, bgi::quadratic<16>>, _NLAYERS_TPC> _rtrees; */
0265   std::array<bgi::rtree<pointKey, bgi::quadratic<16>>, 3> _rtrees;  // need three layers at a time
0266
0267   double Ne_frac = 0.00;
0268   double Ar_frac = 0.75;
0269   double CF4_frac = 0.20;
0270   double N2_frac = 0.00;
0271   double isobutane_frac = 0.05;
0272 };
0273
0274 #endif