packages/trackbase/TrackFitUtils.h

0001 #ifndef TRACKBASE_TRACKFITUTILS_H
0002 #define TRACKBASE_TRACKFITUTILS_H
0003
0004 #include "ActsSurfaceMaps.h"
0005 #include "TrkrDefs.h"
0006 #include "MvtxDefs.h"
0007
0008 #include <Acts/Definitions/Algebra.hpp>
0009
0010 #include <tuple>
0011 #include <utility>
0012 #include <vector>
0013
0014 class ActsGeometry;
0015 class TrkrClusterContainer;
0016
0017 namespace TrackFitUtils
0018 {
0019   using position_t = std::pair<double, double>;
0020   using position_vector_t = std::vector<position_t>;
0021
0022   /// circle fit output [R, x0, y0]
0023   using circle_fit_output_t = std::tuple<double, double, double>;
0024
0025   /**
0026    * Circle fit to a given set of data points (in 2D)
0027    * This is an algebraic fit, due to Taubin, based on the journal article
0028    * G. Taubin, "Estimation Of Planar Curves, Surfaces And Nonplanar
0029    * Space Curves Defined By Implicit Equations, With
0030    * Applications To Edge And Range Image Segmentation",
0031    * IEEE Trans. PAMI, Vol. 13, pages 1115-1138, (1991)
0032    * It works well whether data points are sampled along an entire circle or along a small arc.
0033    * It still has a small bias and its statistical accuracy is slightly lower than that of the geometric fit (minimizing geometric distances),
0034    * It provides a very good initial guess for a subsequent geometric fit.
0035    * Nikolai Chernov  (September 2012)
0036    */
0037   circle_fit_output_t circle_fit_by_taubin(const position_vector_t&);
0038
0039   /// convenient overload
0040   circle_fit_output_t circle_fit_by_taubin(const std::vector<Acts::Vector3>&);
0041
0042   /// line fit output [slope, intercept]
0043   using line_fit_output_t = std::tuple<double, double>;
0044
0045   /// line_fit
0046   /**
0047    * copied from: https://www.bragitoff.com
0048    * typically used to fit we want to fit z vs radius
0049    *
0050    * Updated to use the "Deming model" by default:
0051    * minimizing by orthoncal distance to line in x and y
0052    * (instead of the y-distance). For details, see:
0053    * http://staff.pubhealth.ku.dk/~bxc/MethComp/Deming.pdf
0054    */
0055   line_fit_output_t line_fit(const position_vector_t&);
0056
0057   /*
0058    * Need to make a metric for distance from points to lines origin (pca).
0059    *  - project point "global" to the line.
0060    *  - return distance on line to the pca (the point of closest approach to origin)
0061    */
0062   double line_dist_to_pca (const double slope, const double intercept,
0063       const Acts::Vector2& pca, const Acts::Vector3& global);
0064
0065   /// convenient overload
0066   line_fit_output_t line_fit(const std::vector<Acts::Vector3>&);
0067
0068   line_fit_output_t line_fit_xy(const std::vector<Acts::Vector3>& positions);
0069   line_fit_output_t line_fit_xz(const std::vector<Acts::Vector3>& positions);
0070
0071   /// line-circle intersection output. (xplus, yplus, xminus, yminus)
0072   using line_circle_intersection_output_t = std::tuple<double, double, double, double>;
0073   /**
0074   * r is radius of sPHENIX layer
0075   * m and b are parameters of line fitted to TPC clusters in the x-y plane (slope and intersection)
0076   * the solutions are xplus, xminus, yplus, yminus
0077   * The intersection of the line-circle occurs when
0078   * y = m*x + b
0079   * Here we assume that circle is an sPHENIX layer centered on x1=y1=0
0080   * x^2 +y^2 = r^2
0081   * substitute for y in equation of circle, solve for x, and then for y.
0082   **/
0083   line_circle_intersection_output_t line_circle_intersection(double r, double m, double b);
0084
0085   /// circle-circle intersection output. (xplus, yplus, xminus, yminus)
0086   using circle_circle_intersection_output_t = std::tuple<double, double, double, double>;
0087
0088   /**
0089    * r1 is radius of sPHENIX layer
0090    * r2, x2 and y2 are parameters of circle fitted to TPC clusters
0091    * the solutions are xplus, xminus, yplus, yminus
0092
0093    * The intersection of two circles occurs when
0094    * (x-x1)^2 + (y-y1)^2 = r1^2,  / (x-x2)^2 + (y-y2)^2 = r2^2
0095    * Here we assume that circle 1 is an sPHENIX layer centered on x1=y1=0, and circle 2 is arbitrary
0096    * x^2 +y^2 = r1^2,   (x-x2)^2 + (y-y2)^2 = r2^2
0097    * expand the equations and subtract to eliminate the x^2 and y^2 terms, gives the radical line connecting the intersection points
0098    * iy = - (2*x2*x - D) / 2*y2,
0099    * then substitute for y in equation of circle 1
0100    */
0101   circle_circle_intersection_output_t circle_circle_intersection(double r1, double r2, double x2, double y2);
0102
0103   unsigned int addClusters(std::vector<float>& fitpars,
0104                                   double dca_cut,
0105                                   ActsGeometry* _tGeometry,
0106                                   TrkrClusterContainer* _cluster_map,
0107                                   std::vector<Acts::Vector3>& global_vec,
0108                                   std::vector<TrkrDefs::cluskey>& cluskey_vec,
0109                                   unsigned int startLayer,
0110                                   unsigned int endLayer);
0111
0112   unsigned int addClustersOnLine(TrackFitUtils::line_fit_output_t& fitpars,
0113                                         const bool& isXY,
0114                                         const double& dca_cut,
0115                                         ActsGeometry* tGeometry,
0116                                         TrkrClusterContainer* clusterContainer,
0117                                         std::vector<Acts::Vector3>& global_vec,
0118                                         std::vector<TrkrDefs::cluskey>& cluskey_vec,
0119                                         const unsigned int& startLayer,
0120                                         const unsigned int& endLayer);
0121
0122   std::pair<Acts::Vector3, Acts::Vector3> get_helix_tangent(const std::vector<float>& fitpars, Acts::Vector3& global);
0123
0124   Acts::Vector3 get_helix_pca(std::vector<float>& fitpars, const Acts::Vector3& global);
0125
0126   Acts::Vector2 get_circle_point_pca(float radius, float x0, float y0, Acts::Vector3 global);
0127
0128   std::vector<float> fitClusters(std::vector<Acts::Vector3>& global_vec,
0129                                         const std::vector<TrkrDefs::cluskey> &cluskey_vec,
0130                                         bool use_intt = false);
0131   void getTrackletClusters(ActsGeometry* _tGeometry,
0132                                   TrkrClusterContainer* _cluster_map,
0133                                   std::vector<Acts::Vector3>& global_vec,
0134                                   const std::vector<TrkrDefs::cluskey>& cluskey_vec);
0135   Acts::Vector3 surface_3Dline_intersection(const TrkrDefs::cluskey& key,
0136                                                    TrkrCluster* cluster, ActsGeometry* geometry, float& xyslope, float& xyint, float& yzslope, float& yzint);
0137
0138   Acts::Vector3 getPCALinePoint(const Acts::Vector3& global, const Acts::Vector3& tangent, const Acts::Vector3& posref);
0139   Acts::Vector3 get_helix_surface_intersection(const Surface& surf,
0140                                                       std::vector<float>& fitpars, Acts::Vector3& global, ActsGeometry* _tGeometry);
0141   Acts::Vector3 get_line_plane_intersection(const Acts::Vector3& pca, const Acts::Vector3& tangent,
0142                                                    const Acts::Vector3& sensorCenter, const Acts::Vector3& sensorNormal);
0143
0144   std::vector<double> getLineClusterResiduals(position_vector_t& rz_pts, float slope, float intercept);
0145   std::vector<double> getCircleClusterResiduals(position_vector_t& xy_pts, float R, float X0, float Y0);
0146
0147   Acts::Vector2 get_line_point_pca(double slope, double intercept, Acts::Vector3 global);
0148   std::vector<float> fitClustersZeroField(const std::vector<Acts::Vector3>& global_vec,
0149                                const std::vector<TrkrDefs::cluskey>& cluskey_vec, bool use_intt, bool mvtx_east = false, bool mvtx_west=false);
0150
0151
0152   float get_helix_pathlength(std::vector<float>& fitpars, const Acts::Vector3& start_point, const Acts::Vector3& end_point);
0153   float get_helix_surface_pathlength(const Surface& surf, std::vector<float>& fitpars, const Acts::Vector3& start_point, ActsGeometry* tGeometry);
0154
0155   /* std::tuple<double,double> dca_on_line2D_to_point(const double x0, const double y0, const double xy_m, const double xy_b); */
0156   // get track fit parameters for track matching:: is-good-fit, phi, eta, pt==1, pos_vec3, mom_vec3
0157   std::tuple<bool, double, double, double, Acts::Vector3, Acts::Vector3>
0158     zero_field_track_params(
0159       ActsGeometry* _tGeometry,
0160       TrkrClusterContainer* _cluster_map,
0161       const std::vector<TrkrDefs::cluskey>& clusters
0162     );
0163
0164    double z_fit_to_pca(const double slope, const double intercept,
0165     const std::vector<Acts::Vector3>& glob_pts);
0166   bool isTrackCrossMvtxHalf(const std::vector<TrkrDefs::cluskey> &cluskey_vec);
0167   bool includeMvtxHit(TrkrDefs::cluskey clus_key, bool mvtx_east, bool mvtx_west);
0168   bool isMvtxEast(uint32_t hitsetkey);
0169 };
0170
0171 #endif