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 = default;
0102
0103   void SetSplitSeeds(bool opt = true) { _split_seeds = opt; }
0104   void SetLayerRange(unsigned int layer_low, unsigned int layer_up)
0105   {
0106     _start_layer = layer_low;
0107     _end_layer = layer_up;
0108   }
0109   void SetSearchWindow(float z_width, float phi_width)
0110   {
0111     _neighbor_z_width = z_width;
0112     _neighbor_phi_width = phi_width;
0113   }
0114   void SetClusAdd_delta_window(float _dzdr_cutoff, float _dphidr2_cutoff)
0115   {
0116     _clusadd_delta_dzdr_window = _dzdr_cutoff;
0117     _clusadd_delta_dphidr2_window = _dphidr2_cutoff;
0118   }
0119   void SetMinHitsPerCluster(unsigned int minHits) { _min_nhits_per_cluster = minHits; }
0120   void SetMinClustersPerTrack(unsigned int minClus) { _min_clusters_per_track = minClus; }
0121   void SetNClustersPerSeedRange(unsigned int minClus, unsigned int maxClus)
0122   {
0123     _min_clusters_per_seed = minClus;
0124     _max_clusters_per_seed = maxClus;
0125     if (_min_clusters_per_seed < 3)
0126     {
0127       std::cout << " Error in SetNClustersPerSeedRange: " << __FILE__
0128                 << " min value cannot be less than three." << std::endl;
0129       assert(false);
0130     }
0131   }
0132
0133   // obsolete
0134   void set_field_dir(const double)
0135   { std::cout << "PHCASeeding::set_field_dir - Warning - This function does nothing. Remove from macro" << std::endl; }
0136
0137   // obsolete
0138   void magFieldFile(const std::string&)
0139   { std::cout << "PHCASeeding::magFieldFile - Warning - This function does nothing. Remove from macro" << std::endl; }
0140
0141   // obsolete
0142   void useConstBField(bool)
0143   { std::cout << "PHCASeeding::useConstBField - Warning - This function does nothing. Remove from macro" << std::endl; }
0144
0145   // obsolete
0146   void constBField(float)
0147   { std::cout << "PHCASeeding::constBField - Warning - This function does nothing. Remove from macro" << std::endl; }
0148
0149   void useFixedClusterError(bool opt) { _use_fixed_clus_err = opt; }
0150   void setFixedClusterError(int i, double val) { _fixed_clus_err.at(i) = val; }
0151   void set_pp_mode(bool mode) { _pp_mode = mode; }
0152   void reject_zsize1_clusters(bool mode){_reject_zsize1 = mode;}
0153   void setNeonFraction(double frac) { Ne_frac = frac; };
0154   void setArgonFraction(double frac) { Ar_frac = frac; };
0155   void setCF4Fraction(double frac) { CF4_frac = frac; };
0156   void setNitrogenFraction(double frac) { N2_frac = frac; };
0157   void setIsobutaneFraction(double frac) { isobutane_frac = frac; };
0158
0159  protected:
0160   int Setup(PHCompositeNode* topNode) override;
0161   int Process(PHCompositeNode* topNode) override;
0162   int InitializeGeometry(PHCompositeNode* topNode);
0163   /* int FindSeedsLayerSkip(double cosTheta_limit); */
0164   int End() override;
0165
0166  private:
0167   bool _save_clus_proc = false;
0168   TFile* _f_clustering_process = nullptr;
0169   int _tupout_count = -1;
0170   int _n_tupchains = -1;
0171   // Save the steps of the clustering
0172   TNtuple* _tupclus_all = nullptr;          // all clusters
0173   TNtuple* _tupclus_links = nullptr;        //  clusters which are linked
0174   TNtuple* _tupclus_bilinks = nullptr;      // bi-linked clusters
0175   TNtuple* _tupclus_seeds = nullptr;        // seed (outermost three bi-linked chain
0176   TNtuple* _tupclus_grown_seeds = nullptr;  // seeds saved out
0177                                             //
0178   TNtuple* _tup_chainbody = nullptr;
0179   TNtuple* _tup_chainfork = nullptr;
0180
0181   TNtuple* _tupwin_link = nullptr;       // cluster pairs with dphi and dz (very wide windows) x0,x1,y0,y1,z0,z1
0182   TNtuple* _tupwin_cos_angle = nullptr;  // directed cosine (all cluster triples which are found) x0,x1,x2
0183   TNtuple* _tupwin_seed23 = nullptr;     // from third to fourth link -- can only use passing seeds
0184   TNtuple* _tupwin_seedL1 = nullptr;     // from third to fourth link -- can only use passing seeds
0185
0186   TNtuple* _search_windows = nullptr;  // This is really just a lazy way to store what search paramaters where used.
0187                                        // It would be equally valid in a TMap or map<string,float>
0188   // functions used to fill tuples -- only defined if _PHCASEEDING_CLUSTERLOG_TUPOUT_ is defined in preprocessor
0189   void write_tuples();
0190   void fill_tuple(TNtuple*, float, TrkrDefs::cluskey, const Acts::Vector3&) const;
0191   void fill_tuple_with_seed(TNtuple*, const keyList&, const PositionMap&) const;
0192   void process_tupout_count();
0193   void FillTupWinLink(bgi::rtree<pointKey, bgi::quadratic<16>>&, const coordKey&, const PositionMap&) const;
0194   void FillTupWinCosAngle(const TrkrDefs::cluskey, const TrkrDefs::cluskey, const TrkrDefs::cluskey, const PositionMap&, double cos_angle, bool isneg) const;
0195   void FillTupWinGrowSeed(const keyList& seed, const keyLink& link, const PositionMap& globalPositions) const;
0196   void fill_split_chains(const keyList& chain, const keyList& keylinks, const PositionMap& globalPositions, int& nchains) const;
0197   /* void fill_tuple_with_seed(TN */
0198
0199   // have a comparator to search vector of sorted bilinks for links with starting
0200   // a key of a given value
0201   class CompKeyToBilink
0202   {
0203    public:
0204     bool operator()(const keyLink& p, const TrkrDefs::cluskey& val) const
0205     {
0206       return p.first < val;
0207     }
0208     bool operator()(const TrkrDefs::cluskey& val, const keyLink& p) const
0209     {
0210       return val < p.first;
0211     }
0212   };
0213   std::pair<std::vector<keyLink>::iterator, std::vector<keyLink>::iterator> FindBilinks(const TrkrDefs::cluskey& key);
0214
0215   /// tpc distortion correction utility class
0216   TpcDistortionCorrection m_distortionCorrection;
0217
0218   /// get global position for a given cluster
0219   /**
0220    * uses ActsTransformation to convert cluster local position into global coordinates
0221    * incorporates TPC distortion correction, if present
0222    */
0223   Acts::Vector3 getGlobalPosition(TrkrDefs::cluskey, TrkrCluster*) const;
0224   std::pair<PositionMap, keyListPerLayer> FillGlobalPositions();
0225   std::pair<keyLinks, keyLinkPerLayer> CreateBiLinks(const PositionMap& globalPositions, const keyListPerLayer& ckeys);
0226   PHCASeeding::keyLists FollowBiLinks(const keyLinks& trackSeedPairs, const keyLinkPerLayer& bilinks, const PositionMap& globalPositions) const;
0227   std::vector<coordKey> FillTree(bgi::rtree<pointKey, bgi::quadratic<16>>&, const keyList&, const PositionMap&, int layer);
0228   int FindSeedsWithMerger(const PositionMap&, const keyListPerLayer&);
0229
0230   void QueryTree(const bgi::rtree<pointKey, bgi::quadratic<16>>& rtree, double phimin, double zmin, double phimax, double zmax, std::vector<pointKey>& returned_values) const;
0231   std::vector<TrackSeed_v2> RemoveBadClusters(const std::vector<keyList>& seeds, const PositionMap& globalPositions) const;
0232   double getMengerCurvature(TrkrDefs::cluskey a, TrkrDefs::cluskey b, TrkrDefs::cluskey c, const PositionMap& globalPositions) const;
0233
0234   void publishSeeds(const std::vector<TrackSeed_v2>& seeds) const;
0235
0236   // int _nlayers_all;
0237   // unsigned int _nlayers_seeding;
0238   // std::vector<int> _seeding_layer;
0239
0240   const unsigned int _nlayers_maps;
0241   const unsigned int _nlayers_intt;
0242   const unsigned int _nlayers_tpc;
0243   unsigned int _start_layer;
0244   unsigned int _end_layer;
0245   unsigned int _min_nhits_per_cluster;
0246   unsigned int _min_clusters_per_track;
0247   unsigned int _max_clusters_per_seed = 6;  // currently not used
0248   unsigned int _min_clusters_per_seed = 6;  // currently the abs. number used
0249   //  float _cluster_z_error;
0250   //  float _cluster_alice_y_error;
0251   float _neighbor_phi_width;
0252   float _neighbor_z_width;
0253   float _clusadd_delta_dzdr_window = 0.5;
0254   float _clusadd_delta_dphidr2_window = 0.005;
0255   float _max_sin_phi;
0256   /* float _cosTheta_limit; */
0257   double _rz_outlier_threshold = 0.1;
0258   double _xy_outlier_threshold = 0.1;
0259   bool _split_seeds = true;
0260   bool _reject_zsize1 = false;
0261   bool _use_fixed_clus_err = false;
0262   bool _pp_mode = false;
0263   std::array<double, 3> _fixed_clus_err = {.1, .1, .1};
0264
0265   /// acts geometry
0266   ActsGeometry* m_tGeometry{nullptr};
0267
0268   /// global position wrapper
0269   TpcGlobalPositionWrapper m_globalPositionWrapper;
0270
0271   std::unique_ptr<PHTimer> t_seed;
0272   std::unique_ptr<PHTimer> t_fill;
0273   std::unique_ptr<PHTimer> t_makebilinks;
0274   std::unique_ptr<PHTimer> t_makeseeds;
0275   /* std::array<bgi::rtree<pointKey, bgi::quadratic<16>>, _NLAYERS_TPC> _rtrees; */
0276   std::array<bgi::rtree<pointKey, bgi::quadratic<16>>, 3> _rtrees;  // need three layers at a time
0277
0278   double Ne_frac = 0.00;
0279   double Ar_frac = 0.75;
0280   double CF4_frac = 0.20;
0281   double N2_frac = 0.00;
0282   double isobutane_frac = 0.05;
0283 };
0284
0285 #endif