packages/trackreco/SecondaryVertexFinder.h

0001 // Tell emacs that this is a C++ source
0002 //  -*- C++ -*-.
0003
0004 /*!
0005  *  \file         SecondaryVertexFinder
0006  *  \brief      Class for determining primary vertices usin g fitted tracks
0007  *  \author  Tony Frawley <afrawley@fsu.edu>
0008  */
0009
0010 #ifndef SECONDARYVERTEXFINDER_H
0011 #define SECONDARYVERTEXFINDER_H
0012
0013 #include <fun4all/SubsysReco.h>
0014
0015 #include <trackbase/TpcDefs.h>
0016 #include <trackbase/TrkrDefs.h>
0017
0018 #include <Acts/Definitions/Algebra.hpp>
0019
0020 //#pragma GCC diagnostic push
0021 //#pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0022 #include <Acts/Propagator/Propagator.hpp>
0023 //#pragma GCC diagnostic pop
0024
0025 #include <Acts/EventData/TrackParameters.hpp>
0026 #include <Acts/Surfaces/CylinderSurface.hpp>
0027 #include <Acts/Utilities/Result.hpp>
0028 #include <ActsExamples/EventData/Trajectories.hpp>
0029
0030 #include <map>
0031 #include <set>
0032 #include <string>
0033 #include <vector>
0034
0035 #include <Eigen/Dense>
0036
0037 class PHCompositeNode;
0038 class SvtxTrack;
0039 class SvtxTrackMap;
0040 class SvtxVertexMap;
0041 class TrkrCluster;
0042 class TrackSeed;
0043 class ActsGeometry;
0044 class TrkrClusterContainer;
0045 class TNtuple;
0046 class TH2D;
0047 class TH1D;
0048 class TFile;
0049 class TLorentzVector;
0050
0051 using BoundTrackParam = const Acts::BoundTrackParameters;
0052 using BoundTrackParamResult = Acts::Result<BoundTrackParam>;
0053 using SurfacePtr = std::shared_ptr<const Acts::Surface>;
0054 using Trajectory = ActsExamples::Trajectories;
0055
0056 class SecondaryVertexFinder : public SubsysReco
0057 {
0058  public:
0059   SecondaryVertexFinder(const std::string& name = "SecondaryVertexFinder");
0060
0061   ~SecondaryVertexFinder() override;
0062
0063   int InitRun(PHCompositeNode* topNode) override;
0064   int process_event(PHCompositeNode* topNode) override;
0065   int End(PHCompositeNode* topNode) override;
0066
0067   void setTrackDcaCut(const double cutxy, const double cutz)
0068   {
0069     _track_dcaxy_cut = cutxy;
0070     _track_dcaz_cut = cutz;
0071   }
0072   void setTwoTrackDcaCut(const double cut) { _two_track_dcacut = cut; }
0073   void setMinPathCut(const double cut) { _min_path_cut = cut; }
0074   void setTrackQualityCut(double cut) { _qual_cut = cut; }
0075   void setRequireMVTX(bool set) { _require_mvtx = set; }
0076   void setOutfileName(const std::string& filename) { outfile = filename; }
0077   void setDecayParticleMass(double mass) { _decaymass = mass; }
0078   void set_write_electrons_node(bool flag) { _write_electrons_node = flag; }
0079   void set_write_ntuple(bool flag) { _write_ntuple = flag; }
0080
0081  private:
0082   int GetNodes(PHCompositeNode* topNode);
0083   int CreateOutputNode(PHCompositeNode* topNode);
0084
0085   bool passConversionElectronCuts(const TLorentzVector& tsum,
0086                                   SvtxTrack* tr1, SvtxTrack* tr2, float pair_dca,
0087                                   const Eigen::Vector3d& PCA, const Eigen::Vector3d& VTX);
0088
0089   bool hasSiliconSeed(SvtxTrack* tr);
0090   void outputTrackDetails(SvtxTrack* tr);
0091   void get_dca(SvtxTrack* track, float& dca3dxy, float& dca3dz, float& dca3dxysigma, float& dca3dzsigma);
0092   bool circle_circle_intersection(double r0, double x0, double y0, double r1, double x1, double y1, std::vector<double>& intersectionXY);
0093
0094   void findPcaTwoLines(Eigen::Vector3d pos1, Eigen::Vector3d mom1, Eigen::Vector3d pos2, Eigen::Vector3d mom2,
0095                        double& dca, Eigen::Vector3d& PCA1, Eigen::Vector3d& PCA2);
0096   void getCircleXYTrack(SvtxTrack* track, double& R, Eigen::Vector2d& center);
0097   double getZFromIntersectionXY(SvtxTrack* track, double& R, Eigen::Vector2d& center, Eigen::Vector2d intersection);
0098   bool projectTrackToPoint(SvtxTrack* track, Eigen::Vector3d& PCA, Eigen::Vector3d& pos, Eigen::Vector3d& mom);
0099
0100   void fillNtp(SvtxTrack* track1, SvtxTrack* track2, double dca3dxy1, double dca3dz1, double dca3dxy2, double dca3dz2, Eigen::Vector3d vpos1, Eigen::Vector3d vmom1, Eigen::Vector3d vpos2, Eigen::Vector3d vmom2, Acts::Vector3 pca_rel1, Acts::Vector3 pca_rel2, double pair_dca, double invariantMass, double invariantPt, double path, int has_silicon_1, int has_siilicon_2);
0101
0102   SvtxTrackMap* _track_map{nullptr};
0103   SvtxTrackMap* _track_map_electrons{nullptr};
0104   //  SvtxTrack *_track{nullptr};
0105   SvtxVertexMap* _svtx_vertex_map{nullptr};
0106   ActsGeometry* _tGeometry{nullptr};
0107
0108   bool _require_mvtx = false;
0109   bool _write_electrons_node = true;
0110   bool _write_ntuple = false;
0111
0112   double _decaymass = 0.000511;  // conversion electrons, default
0113
0114   // these are minimal cuts used to make the ntuple
0115   // They can be tightened later when analyzing the ntuple
0116
0117   // single track cuts
0118   double _track_dcaxy_cut = 0.020;
0119   double _track_dcaz_cut = 0.020;
0120   double _qual_cut = 4.0;
0121
0122   // track_pair cuts
0123   double _two_track_dcacut = 0.5;          // 5000 microns
0124   double _max_intersection_radius = 40.0;  // discard intersections at greater than 40 cm radius
0125   double _projected_track_z_cut = 1.0;
0126
0127   // decay vertex cuts
0128   double _min_path_cut = 0.2;
0129   double _costheta_cut = 0.9985;
0130
0131   // specific conversion electron cuts
0132   double _conversion_pair_dcacut = 0.2;  // 2000 microns
0133   unsigned int _min_tpc_clusters = 40;
0134   double _deta_cut = 0.05;
0135   double _invariant_pt_cut = 0.1;
0136   double _max_mass_cut = 0.03;
0137
0138   TH2D* recomass{nullptr};
0139   TH2D* hdecaypos{nullptr};
0140   TH1D* hdecay_radius{nullptr};
0141   TNtuple* ntp{nullptr};
0142   std::string outfile;
0143 };
0144
0145 #endif  // SECONDARYVERTEXFINDER_H