packages/trackreco/PHSimpleVertexFinder.h

0001 // Tell emacs that this is a C++ source
0002 //  -*- C++ -*-.
0003
0004 /*!
0005  *  \file         PHSimpleVertexFinder
0006  *  \brief      Class for determining primary vertices usin g fitted tracks
0007  *  \author  Tony Frawley <afrawley@fsu.edu>
0008  */
0009
0010 #ifndef PHSIMPLEVERTEXFINDER_H
0011 #define PHSIMPLEVERTEXFINDER_H
0012
0013 #include <fun4all/SubsysReco.h>
0014 #include <trackbase/TrkrDefs.h>
0015 #include <trackbase/ActsGeometry.h>
0016
0017 #include <map>
0018 #include <set>
0019 #include <string>
0020 #include <vector>
0021
0022 #include <Eigen/Dense>
0023
0024 class PHCompositeNode;
0025 class SvtxTrack;
0026 class SvtxTrackMap;
0027 class SvtxVertexMap;
0028 class TrkrCluster;
0029 class TrackVertexCrossingAssoc;
0030 class TrackSeed;
0031 class TrkrClusterContainer;
0032 class ActsGeometry;
0033
0034 class PHSimpleVertexFinder : public SubsysReco
0035 {
0036  public:
0037   PHSimpleVertexFinder(const std::string &name = "PHSimpleVertexFinder");
0038
0039   ~PHSimpleVertexFinder() override = default;
0040
0041   int InitRun(PHCompositeNode *topNode) override;
0042   int process_event(PHCompositeNode *topNode) override;
0043   int End(PHCompositeNode *topNode) override;
0044   void setBeamLineCut(const double cut) { _beamline_xy_cut = cut; }
0045   void setDcaCut(const double cut) { _base_dcacut = cut; }
0046   void setTrackQualityCut(double cut) { _qual_cut = cut; }
0047   void setRequireMVTX(bool set) { _require_mvtx = set; }
0048   void setNmvtxRequired(unsigned int n) { _nmvtx_required = n; }
0049   void setTrackPtCut(const double cut) { _track_pt_cut = cut; }
0050   // void setUseTrackCovariance(bool set) {_use_track_covariance = set;}
0051   void setOutlierPairCut(const double cut) { _outlier_cut = cut; }
0052   void setTrackMapName(const std::string &name) { _track_map_name = name; }
0053   void setVertexMapName(const std::string &name) { _vertex_map_name = name; }
0054   void zeroField(const bool flag) { _zero_field = flag; }
0055   void setTrkrClusterContainerName(std::string &name){ m_clusterContainerName = name; }
0056   void set_pp_mode(bool mode) { _pp_mode = mode; }
0057
0058  private:
0059   int GetNodes(PHCompositeNode *topNode);
0060   int CreateNodes(PHCompositeNode *topNode);
0061
0062   void checkDCAs(SvtxTrackMap *track_map);
0063   void checkDCAsZF(SvtxTrackMap *track_map);
0064   void checkDCAs();
0065
0066   void getTrackletClusterList(TrackSeed* tracklet, std::vector<TrkrDefs::cluskey>& cluskey_vec);
0067
0068   void findDcaTwoTracks(SvtxTrack *tr1, SvtxTrack *tr2);
0069   double dcaTwoLines(const Eigen::Vector3d &p1, const Eigen::Vector3d &v1,
0070                      const Eigen::Vector3d &p2, const Eigen::Vector3d &v2,
0071                      Eigen::Vector3d &PCA1, Eigen::Vector3d &PCA2);
0072   std::vector<std::set<unsigned int>> findConnectedTracks();
0073   void removeOutlierTrackPairs();
0074   double getMedian(std::vector<double> &v);
0075   double getAverage(std::vector<double> &v);
0076
0077   SvtxTrackMap *_track_map{nullptr};
0078   TrkrClusterContainer* _cluster_map{nullptr};
0079   SvtxVertexMap *_svtx_vertex_map{nullptr};
0080   ActsGeometry* _tGeometry{nullptr};
0081
0082   double _base_dcacut = 0.0080;  // 80 microns
0083   double _active_dcacut = 0.080;
0084   double _beamline_xy_cut = 0.2;  // must be within 2 mm of beam line
0085   double _qual_cut = 10.0;
0086   bool _require_mvtx = true;
0087   unsigned int _nmvtx_required = 3;
0088   double _track_pt_cut = 0.0;
0089   double _outlier_cut = 0.015;
0090   //name of TRKR_CLUSTER Container
0091   std::string m_clusterContainerName = "TRKR_CLUSTER";
0092
0093   bool _zero_field = false;     // fit straight lines if true
0094
0095   std::string _track_map_name = "SvtxTrackMap";
0096   std::string _vertex_map_name = "SvtxVertexMap";
0097   std::multimap<unsigned int, unsigned int> _vertex_track_map;
0098   using matrix_t = Eigen::Matrix<double, 3, 3>;
0099   std::multimap<unsigned int, std::pair<unsigned int, double>> _track_pair_map;
0100   // Eigen::Vector3d is an Eigen::Matrix<double,3,1>
0101   std::multimap<unsigned int, std::pair<unsigned int, std::pair<Eigen::Vector3d,
0102                                                                 Eigen::Vector3d>>>
0103       _track_pair_pca_map;
0104   std::map<unsigned int, Eigen::Vector3d> _vertex_position_map;
0105   std::map<unsigned int, matrix_t> _vertex_covariance_map;
0106   std::set<unsigned int> _vertex_set;
0107
0108   TrackVertexCrossingAssoc *_track_vertex_crossing_map{nullptr};
0109
0110   bool _pp_mode = true;  // default to pp mode
0111 };
0112
0113 #endif  // PHSIMPLEVERTEXFINDER_H