packages/trackbase_historic/WeightedTrackZeroField.cc

0001 #include "WeightedTrackZeroField.h"
0002
0003 #include <phool/phool.h>
0004
0005 #include <stdexcept>
0006
0007 void
0008 WeightedTrackZeroField::set_parameters (
0009     double const* params
0010 ) {
0011     m_cos_theta = std::cos(params[2]);
0012     m_sin_theta = std::sin(params[2]);
0013     m_cos_phi = std::cos(params[3]);
0014     m_sin_phi = std::sin(params[3]);
0015
0016     m_intercept = Eigen::Vector3d {
0017         params[0], // x-component of xy plane intercept
0018         params[1], // y-component of xy plane intercept
0019         0,
0020     };
0021
0022     m_slope = Eigen::Vector3d {
0023         m_sin_theta * m_cos_phi,
0024         m_sin_theta * m_sin_phi,
0025         m_cos_theta
0026     };
0027 }
0028
0029 void
0030 WeightedTrackZeroField::get_parameters (
0031     double* params
0032 ) const {
0033     params[0] = m_intercept(0);
0034     params[1] = m_intercept(1);
0035     params[2] = std::acos(m_slope(2));
0036     params[3] = std::atan2(m_slope(1), m_slope(0));
0037 }
0038
0039 void
0040 WeightedTrackZeroField::configure_minimizer (
0041     ROOT::Math::Minimizer& minimizer
0042 ) {
0043     if (size() < 2) {
0044         std::stringstream what;
0045         what
0046             << PHWHERE
0047             << " too few points to perform fit"
0048             << std::endl;
0049         throw std::runtime_error(what.str());
0050     }
0051
0052     // Initial slope
0053     m_slope = (operator[](1).cluster_position - operator[](0).cluster_position).normalized();
0054     double phi = std::atan2(m_slope(1), m_slope(0));
0055     double theta = std::acos(m_slope(2));
0056     m_cos_theta = m_slope(2);
0057     m_sin_theta = std::sqrt(m_slope(0) * m_slope(0) + m_slope(1) * m_slope(1));
0058     m_cos_phi = m_slope(0) / m_sin_theta;
0059     m_sin_phi = m_slope(1) / m_sin_theta;
0060
0061     // Initial intercept
0062     // Solve for s:
0063     // 0 = (m_slope * s + operator[](0).cluster_position).dot({0, 0, 1})
0064     double s = -operator[](0).cluster_position(2) / m_slope(2);
0065     m_intercept = m_slope * s + operator[](0).cluster_position;
0066
0067     minimizer.SetVariable(0, "x",     m_intercept(0), 1E-6);
0068     minimizer.SetVariable(1, "y",     m_intercept(1), 1E-6);
0069     minimizer.SetVariable(2, "theta", theta,          1E-6);
0070     minimizer.SetVariable(3, "phi",   phi,            1E-6);
0071
0072     // Bounds for angular parameters
0073     minimizer.SetVariableLimits(2,          0,     3.1416);
0074     minimizer.SetVariableLimits(3, phi - 1.58, phi + 1.58);
0075 }
0076
0077 Eigen::Vector3d
0078 WeightedTrackZeroField::get_partial_derivative (
0079     int index,
0080     double path_length
0081 ) const {
0082     switch (index) {
0083         case 0: return Eigen::Vector3d {1.0, 0.0, 0.0};
0084         case 1: return Eigen::Vector3d {0.0, 1.0, 0.0};
0085         case 2: return Eigen::Vector3d {  m_cos_theta * m_cos_phi, m_cos_theta * m_sin_phi, -m_sin_theta} * path_length;
0086         case 3: return Eigen::Vector3d { -m_sin_theta * m_sin_phi, m_sin_theta * m_cos_phi, 0.0} * path_length;
0087         default: return Eigen::Vector3d::Zero();
0088     }
0089 }
0090
0091 Eigen::Vector3d
0092 WeightedTrackZeroField::get_position_at_path_length (
0093     double path_length
0094 ) const {
0095     return m_slope * path_length + m_intercept;
0096 }
0097
0098 Eigen::Vector3d
0099 WeightedTrackZeroField::get_slope_at_path_length (
0100     double /*path_length*/
0101 ) const {
0102     return m_slope;
0103 }
0104
0105 double
0106 WeightedTrackZeroField::get_path_length_of_intersection (
0107     Eigen::Affine3d const& plane
0108 ) const {
0109     // Solve s:
0110     // 0 = (m_slope * s + m_intercept - plane.translation()).dot(plane.rotation().col(2))
0111     return (plane.translation() - m_intercept).dot(plane.rotation().col(2)) / m_slope.dot(plane.rotation().col(2));
0112 }
0113
0114 double
0115 WeightedTrackZeroField::get_path_length_of_pca (
0116     Eigen::Vector3d const& point
0117 ) const {
0118     // Solve s:
0119     // 0 = (m_slope * s + m_intercept - point).dot(m_slope)
0120     return (point - m_intercept).dot(m_slope);
0121 }
0122