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 setBeamSpotCutX(const double cutlo, const double cuthi) { _beamline_x_cut_lo = cutlo; _beamline_x_cut_hi = cuthi; }
0046   void setBeamSpotCutY(const double cutlo, const double cuthi) { _beamline_y_cut_lo = cutlo; _beamline_y_cut_hi = cuthi; }
0047   void setDcaCut(const double cut) { _base_dcacut = cut; }
0048   void setTrackQualityCut(double cut) { _qual_cut = cut; }
0049   void setRequireMVTX(bool set) { _require_mvtx = set; }
0050   void setNmvtxRequired(unsigned int n) { _nmvtx_required = n; }
0051   void setTrackPtCut(const double cut) { _track_pt_cut = cut; }
0052   // void setUseTrackCovariance(bool set) {_use_track_covariance = set;}
0053   void setOutlierPairCut(const double cut) { _outlier_cut = cut; }
0054   void setTrackMapName(const std::string &name) { _track_map_name = name; }
0055   void setVertexMapName(const std::string &name) { _vertex_map_name = name; }
0056   void zeroField(const bool flag) { _zero_field = flag; }
0057   void setTrkrClusterContainerName(std::string &name){ m_clusterContainerName = name; }
0058   void set_pp_mode(bool mode) { _pp_mode = mode; }
0059
0060  private:
0061   int GetNodes(PHCompositeNode *topNode);
0062   int CreateNodes(PHCompositeNode *topNode);
0063
0064   void checkDCAs(SvtxTrackMap *track_map);
0065   void checkDCAsZF(SvtxTrackMap *track_map);
0066   void checkDCAs();
0067
0068   void getTrackletClusterList(TrackSeed* tracklet, std::vector<TrkrDefs::cluskey>& cluskey_vec);
0069
0070   void findDcaTwoTracks(SvtxTrack *tr1, SvtxTrack *tr2);
0071   double dcaTwoLines(const Eigen::Vector3d &a1, const Eigen::Vector3d &b1,
0072                      const Eigen::Vector3d &a2, const Eigen::Vector3d &b2,
0073                      Eigen::Vector3d &PCA1, Eigen::Vector3d &PCA2);
0074   std::vector<std::set<unsigned int>> findConnectedTracks();
0075   void removeOutlierTrackPairs();
0076   double getMedian(std::vector<double> &v);
0077   double getAverage(std::vector<double> &v);
0078
0079   SvtxTrackMap *_track_map{nullptr};
0080   TrkrClusterContainer* _cluster_map{nullptr};
0081   SvtxVertexMap *_svtx_vertex_map{nullptr};
0082   ActsGeometry* _tGeometry{nullptr};
0083
0084   double _base_dcacut = 0.05;  // pair dca cut - 1000 microns
0085   double _active_dcacut = 0.05;
0086   double _beamline_xy_cut = 0.2;  // must be within this distance of beam line - no longer used
0087   // defines a box around the beam spot
0088   double _beamline_x_cut_lo = -0.2;
0089   double _beamline_x_cut_hi = 0.2;
0090   double _beamline_y_cut_lo = -0.2;
0091   double _beamline_y_cut_hi = 0.2;
0092   double _qual_cut = 10.0;
0093   bool _require_mvtx = true;
0094   unsigned int _nmvtx_required = 3;
0095   double _track_pt_cut = 0.0;
0096   double _outlier_cut = 0.015;
0097
0098   //name of TRKR_CLUSTER Container
0099   std::string m_clusterContainerName = "TRKR_CLUSTER";
0100
0101   bool _zero_field = false;     // fit straight lines if true
0102
0103   std::string _track_map_name = "SvtxTrackMap";
0104   std::string _vertex_map_name = "SvtxVertexMap";
0105   std::multimap<unsigned int, unsigned int> _vertex_track_map;
0106   using matrix_t = Eigen::Matrix<double, 3, 3>;
0107   std::multimap<unsigned int, std::pair<unsigned int, double>> _track_pair_map;
0108   // Eigen::Vector3d is an Eigen::Matrix<double,3,1>
0109   std::multimap<unsigned int, std::pair<unsigned int, std::pair<Eigen::Vector3d,
0110                                                                 Eigen::Vector3d>>>
0111       _track_pair_pca_map;
0112   std::map<unsigned int, Eigen::Vector3d> _vertex_position_map;
0113   std::map<unsigned int, matrix_t> _vertex_covariance_map;
0114   std::set<unsigned int> _vertex_set;
0115
0116   TrackVertexCrossingAssoc *_track_vertex_crossing_map{nullptr};
0117
0118   bool _pp_mode = true;  // default to pp mode
0119 };
0120
0121 #endif  // PHSIMPLEVERTEXFINDER_H