Skip to content

Commit

Permalink
Rework the Riemann fit and broken line fit (#338)
Browse files Browse the repository at this point in the history 
Merge the Riemann and broken line fits into single configurable pixel
n-tuplet fitter, and extend it to work with up to 5 hits.

Mmake the broken line fit the default algorithm.

Try both triplets and quadruplets in the pixel "hole".

Limit pT used to compute the multple scattering.

Use the inline Cholesky decomposition.

Generic clean up and improvements.
  • Loading branch information
@VinInn @fwyzard
VinInn authored and fwyzard committed Jan 15, 2021
1 parent f838bac commit 09a3ecc
Show file tree
Hide file tree
Showing 24 changed files with 1,784 additions and 1,030 deletions.
1,144 changes: 539 additions & 605 deletions RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h
View file

Large diffs are not rendered by default.

93 changes: 27 additions & 66 deletions RecoPixelVertexing/PixelTrackFitting/interface/FitResult.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,97 +2,58 @@
#define RecoPixelVertexing_PixelTrackFitting_interface_FitResult_h

#include <cmath>
#include <cstdint>

#include <cuda_runtime.h>
#include <Eigen/Core>
#include <Eigen/Eigenvalues>

#include "HeterogeneousCore/CUDAUtilities/interface/cuda_assert.h"

namespace Rfit
{

constexpr double d = 1.e-4; //!< used in numerical derivative (J2 in Circle_fit())
constexpr unsigned int max_nop = 4; //!< In order to avoid use of dynamic memory

using Vector2d = Eigen::Vector2d;
using Vector3d = Eigen::Vector3d;
using Vector4d = Eigen::Vector4d;
using Vector5d = Eigen::Matrix<double, 5, 1>;
using Matrix2d = Eigen::Matrix2d;
using Matrix3d = Eigen::Matrix3d;
using Matrix4d = Eigen::Matrix4d;
using Matrix5d = Eigen::Matrix<double, 5, 5>;
using Matrix6d = Eigen::Matrix<double, 6, 6>;

using VectorXd = Eigen::VectorXd;
using MatrixXd = Eigen::MatrixXd;
template<int N>
using MatrixNd = Eigen::Matrix<double, N, N>;
template<int N>
using ArrayNd = Eigen::Array<double, N, N>;
template<int N>
using Matrix2Nd = Eigen::Matrix<double, 2 * N, 2 * N>;
template<int N>
using Matrix3Nd = Eigen::Matrix<double, 3 * N, 3 * N>;
template<int N>
using Matrix2xNd = Eigen::Matrix<double, 2, N>;
template<int N>
using Array2xNd = Eigen::Array<double, 2, N>;
template<int N>
using Matrix3xNd = Eigen::Matrix<double, 3, N>;
template<int N>
using MatrixNx3d = Eigen::Matrix<double, N, 3>;
template<int N>
using MatrixNx5d = Eigen::Matrix<double, N, 5>;
template<int N>
using VectorNd = Eigen::Matrix<double, N, 1>;
template<int N>
using Vector2Nd = Eigen::Matrix<double, 2 * N, 1>;
template<int N>
using Vector3Nd = Eigen::Matrix<double, 3 * N, 1>;
template<int N>
using RowVectorNd = Eigen::Matrix<double, 1, 1, N>;
template<int N>
using RowVector2Nd = Eigen::Matrix<double, 1, 2 * N>;

template<int N>
using Matrix3xNd = Eigen::Matrix<double, 3, N>; // used for inputs hits


using Vector2d = Eigen::Vector2d;
using Vector3d = Eigen::Vector3d;
using Vector4d = Eigen::Vector4d;
using Matrix2d = Eigen::Matrix2d;
using Matrix3d = Eigen::Matrix3d;
using Matrix4d = Eigen::Matrix4d;
using Matrix5d = Eigen::Matrix<double, 5, 5>;
using Matrix6d = Eigen::Matrix<double, 6, 6>;
using Vector5d = Eigen::Matrix<double, 5, 1>;

using Matrix3f = Eigen::Matrix3f;
using Vector3f = Eigen::Vector3f;
using Vector4f = Eigen::Vector4f;
using Vector6f = Eigen::Matrix<double, 6, 1>;

using u_int = unsigned int;


struct circle_fit
{
struct circle_fit
{
Vector3d par; //!< parameter: (X0,Y0,R)
Matrix3d cov;
/*!< covariance matrix: \n
|cov(X0,X0)|cov(Y0,X0)|cov( R,X0)| \n
|cov(X0,Y0)|cov(Y0,Y0)|cov( R,Y0)| \n
|cov(X0, R)|cov(Y0, R)|cov( R, R)|
*/
*/
int32_t q; //!< particle charge
float chi2 = 0.0;
};

struct line_fit
{
};
struct line_fit
{
Vector2d par; //!<(cotan(theta),Zip)
Matrix2d cov;
/*!<
|cov(c_t,c_t)|cov(Zip,c_t)| \n
|cov(c_t,Zip)|cov(Zip,Zip)|
*/
*/
double chi2 = 0.0;
};

struct helix_fit
{
};
struct helix_fit
{
Vector5d par; //!<(phi,Tip,pt,cotan(theta)),Zip)
Matrix5d cov;
/*!< ()->cov() \n
Expand All @@ -101,12 +62,12 @@ struct helix_fit
|(phi,p_t)|(Tip,p_t)|(p_t,p_t)|(c_t,p_t)|(Zip,p_t)| \n
|(phi,c_t)|(Tip,c_t)|(p_t,c_t)|(c_t,c_t)|(Zip,c_t)| \n
|(phi,Zip)|(Tip,Zip)|(p_t,Zip)|(c_t,Zip)|(Zip,Zip)|
*/
*/
float chi2_circle;
float chi2_line;
// Vector4d fast_fit;
// Vector4d fast_fit;
int32_t q; //!< particle charge
} __attribute__((aligned(16)));
}; // __attribute__((aligned(16)));

} // namespace RFit
#endif
203 changes: 203 additions & 0 deletions RecoPixelVertexing/PixelTrackFitting/interface/FitUtils.h
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,203 @@
#ifndef RecoPixelVertexing_PixelTrackFitting_interface_FitUtils_h
#define RecoPixelVertexing_PixelTrackFitting_interface_FitUtils_h


#include "RecoPixelVertexing/PixelTrackFitting/interface/FitResult.h"

#include "HeterogeneousCore/CUDAUtilities/interface/cuda_assert.h"

#include "DataFormats/Math/interface/choleskyInversion.h"


namespace Rfit
{


constexpr double d = 1.e-4; //!< used in numerical derivative (J2 in Circle_fit())



using VectorXd = Eigen::VectorXd;
using MatrixXd = Eigen::MatrixXd;
template<int N>
using MatrixNd = Eigen::Matrix<double, N, N>;
template<int N>
using MatrixNplusONEd = Eigen::Matrix<double, N+1, N+1>;
template<int N>
using ArrayNd = Eigen::Array<double, N, N>;
template<int N>
using Matrix2Nd = Eigen::Matrix<double, 2 * N, 2 * N>;
template<int N>
using Matrix3Nd = Eigen::Matrix<double, 3 * N, 3 * N>;
template<int N>
using Matrix2xNd = Eigen::Matrix<double, 2, N>;
template<int N>
using Array2xNd = Eigen::Array<double, 2, N>;
template<int N>
using MatrixNx3d = Eigen::Matrix<double, N, 3>;
template<int N>
using MatrixNx5d = Eigen::Matrix<double, N, 5>;
template<int N>
using VectorNd = Eigen::Matrix<double, N, 1>;
template<int N>
using VectorNplusONEd = Eigen::Matrix<double, N+1, 1>;
template<int N>
using Vector2Nd = Eigen::Matrix<double, 2 * N, 1>;
template<int N>
using Vector3Nd = Eigen::Matrix<double, 3 * N, 1>;
template<int N>
using RowVectorNd = Eigen::Matrix<double, 1, 1, N>;
template<int N>
using RowVector2Nd = Eigen::Matrix<double, 1, 2 * N>;


using Matrix2x3d = Eigen::Matrix<double, 2, 3>;


using Matrix3f = Eigen::Matrix3f;
using Vector3f = Eigen::Vector3f;
using Vector4f = Eigen::Vector4f;
using Vector6f = Eigen::Matrix<double, 6, 1>;




using u_int = unsigned int;




template <class C>
__host__ __device__ void printIt(C* m, const char* prefix = "")
{
#ifdef RFIT_DEBUG
for (u_int r = 0; r < m->rows(); ++r)
{
for (u_int c = 0; c < m->cols(); ++c)
{
printf("%s Matrix(%d,%d) = %g\n", prefix, r, c, (*m)(r, c));
}
}
#endif
}

/*!
\brief raise to square.
*/
template <typename T>
constexpr T sqr(const T a)
{
return a * a;
}

/*!
\brief Compute cross product of two 2D vector (assuming z component 0),
returning z component of the result.
\param a first 2D vector in the product.
\param b second 2D vector in the product.
\return z component of the cross product.
*/

__host__ __device__ inline double cross2D(const Vector2d& a, const Vector2d& b)
{
return a.x() * b.y() - a.y() * b.x();
}

/*!
* load error in CMSSW format to our formalism
*
*/
template<typename M6xNf, typename M2Nd>
__host__ __device__ void loadCovariance2D(M6xNf const & ge, M2Nd & hits_cov) {
// Index numerology:
// i: index of the hits/point (0,..,3)
// j: index of space component (x,y,z)
// l: index of space components (x,y,z)
// ge is always in sync with the index i and is formatted as:
// ge[] ==> [xx, xy, yy, xz, yz, zz]
// in (j,l) notation, we have:
// ge[] ==> [(0,0), (0,1), (1,1), (0,2), (1,2), (2,2)]
// so the index ge_idx corresponds to the matrix elements:
// | 0 1 3 |
// | 1 2 4 |
// | 3 4 5 |
constexpr uint32_t hits_in_fit = M6xNf::ColsAtCompileTime;
for (uint32_t i=0; i< hits_in_fit; ++i) {
auto ge_idx = 0; auto j=0; auto l=0;
hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) = ge.col(i)[ge_idx];
ge_idx = 2; j=1; l=1;
hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) = ge.col(i)[ge_idx];
ge_idx = 1; j=1; l=0;
hits_cov(i + l * hits_in_fit, i + j * hits_in_fit) =
hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) = ge.col(i)[ge_idx];
}
}

template<typename M6xNf, typename M3xNd>
__host__ __device__ void loadCovariance(M6xNf const & ge, M3xNd & hits_cov) {

// Index numerology:
// i: index of the hits/point (0,..,3)
// j: index of space component (x,y,z)
// l: index of space components (x,y,z)
// ge is always in sync with the index i and is formatted as:
// ge[] ==> [xx, xy, yy, xz, yz, zz]
// in (j,l) notation, we have:
// ge[] ==> [(0,0), (0,1), (1,1), (0,2), (1,2), (2,2)]
// so the index ge_idx corresponds to the matrix elements:
// | 0 1 3 |
// | 1 2 4 |
// | 3 4 5 |
constexpr uint32_t hits_in_fit = M6xNf::ColsAtCompileTime;
for (uint32_t i=0; i<hits_in_fit; ++i) {
auto ge_idx = 0; auto j=0; auto l=0;
hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) = ge.col(i)[ge_idx];
ge_idx = 2; j=1; l=1;
hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) = ge.col(i)[ge_idx];
ge_idx = 5; j=2; l=2;
hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) = ge.col(i)[ge_idx];
ge_idx = 1; j=1; l=0;
hits_cov(i + l * hits_in_fit, i + j * hits_in_fit) =
hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) = ge.col(i)[ge_idx];
ge_idx = 3; j=2; l=0;
hits_cov(i + l * hits_in_fit, i + j * hits_in_fit) =
hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) = ge.col(i)[ge_idx];
ge_idx = 4; j=2; l=1;
hits_cov(i + l * hits_in_fit, i + j * hits_in_fit) =
hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) = ge.col(i)[ge_idx];
}
}

/*!
\brief Transform circle parameter from (X0,Y0,R) to (phi,Tip,p_t) and
consequently covariance matrix.
\param circle_uvr parameter (X0,Y0,R), covariance matrix to
be transformed and particle charge.
\param B magnetic field in Gev/cm/c unit.
\param error flag for errors computation.
*/
__host__ __device__
inline void par_uvrtopak(circle_fit& circle, const double B, const bool error)
{
Vector3d par_pak;
const double temp0 = circle.par.head(2).squaredNorm();
const double temp1 = sqrt(temp0);
par_pak << atan2(circle.q * circle.par(0), -circle.q * circle.par(1)),
circle.q * (temp1 - circle.par(2)), circle.par(2) * B;
if (error)
{
const double temp2 = sqr(circle.par(0)) * 1. / temp0;
const double temp3 = 1. / temp1 * circle.q;
Matrix3d J4;
J4 << -circle.par(1) * temp2 * 1. / sqr(circle.par(0)), temp2 * 1. / circle.par(0), 0.,
circle.par(0) * temp3, circle.par(1) * temp3, -circle.q,
0., 0., B;
circle.cov = J4 * circle.cov * J4.transpose();
}
circle.par = par_pak;
}

}


#endif
27 changes: 27 additions & 0 deletions RecoPixelVertexing/PixelTrackFitting/interface/PixelNtupletsFitter.h
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@
#ifndef RecoPixelVertexing_PixelTrackFitting_PixelNtupletsFitter_H
#define RecoPixelVertexing_PixelTrackFitting_PixelNtupletsFitter_H

#include <vector>

#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "RecoPixelVertexing/PixelTrackFitting/interface/PixelFitterBase.h"
#include "RecoTracker/TkTrackingRegions/interface/TrackingRegion.h"

class PixelNtupletsFitter final : public PixelFitterBase {
public:
explicit PixelNtupletsFitter(float nominalB, const MagneticField *field, bool useRiemannFit);
~PixelNtupletsFitter() override = default;
std::unique_ptr<reco::Track> run(const std::vector<const TrackingRecHit *>& hits,
const TrackingRegion& region,
const edm::EventSetup& setup) const override;

private:
float nominalB_;
const MagneticField *field_;
bool useRiemannFit_;
};

#endif // RecoPixelVertexing_PixelTrackFitting_PixelNtupletsFitter_H
Loading

0 comments on commit 09a3ecc

Please sign in to comment.