FastCircleFit Class Reference

#include <FastCircleFit.h>

Public Member Functions
float	chi2 () const
template<typename P , typename E >
	FastCircleFit (const P &points, const E &errors)
template<size_t N>
	FastCircleFit (const std::array< GlobalPoint, N > &points, const std::array< GlobalError, N > &errors)
float	rho () const
float	x0 () const
float	y0 () const
	~FastCircleFit ()=default

Private Member Functions
template<typename P , typename E , typename C >
void	calculate (const P &points, const E &errors, C &x, C &y, C &z, C &weight)
template<typename T >
T	sqr (T t)

Private Attributes
float	chi2_
float	rho_
float	x0_
float	y0_

Detailed Description

Same (almost) as FastCircle but with arbitrary number of hits.

Strandlie, Wroldsen, Frühwirth NIM A 488 (2002) 332-341. Frühwirth, Strandlie, Waltenberger NIM A 490 (2002) 366-378.

Definition at line 27 of file FastCircleFit.h.

Constructor & Destructor Documentation

template<typename P , typename E >

FastCircleFit::FastCircleFit ( const P &  points, const E &  errors  )

inline

Constructor for containers of GlobalPoint and GlobalError

Template Parameters

P	Container of GlobalPoint
E	Container of GlobalError

Container can be e.g. std::vector or DynArray.

Definition at line 38 of file FastCircleFit.h.

References calculate(), declareDynArray, N, histoStyle::weight, x, y, and z.

                                                  {
    const auto N = points.size();
    declareDynArray(float, N, x);
    declareDynArray(float, N, y);
    declareDynArray(float, N, z);
    declareDynArray(float, N, weight);  // 1/sigma^2
    calculate(points, errors, x, y, z, weight);
  }

template<size_t N>

FastCircleFit::FastCircleFit ( const std::array< GlobalPoint, N > &  points, const std::array< GlobalError, N > &  errors  )

inline

Constructor for std::array of GlobalPoint and GlobalError

Definition at line 51 of file FastCircleFit.h.

References calculate(), histoStyle::weight, x, y, and z.

                                                                                                {
    std::array<float, N> x;
    std::array<float, N> y;
    std::array<float, N> z;
    std::array<float, N> weight;  // 1/sigma^2
    calculate(points, errors, x, y, z, weight);
  }

FastCircleFit::~FastCircleFit ( )

default

Member Function Documentation

template<typename P , typename E , typename C >

void FastCircleFit::calculate ( const P &  points, const E &  errors, C &  x, C &  y, C &  z, C &  weight  )

private

Definition at line 84 of file FastCircleFit.h.

References funct::A, c, chi2_, HLT_FULL_cff::distance, mps_fire::i, SiStripPI::mean, N, dqmiodumpmetadata::n, AlCaHLTBitMon_ParallelJobs::p, point, PRINT, rho_, S(), sqr(), mathSSE::sqrt(), createJobs::tmp, w, x0_, and y0_.

Referenced by FastCircleFit().

                                                                                           {
  const auto N = points.size();

  // transform
  for (size_t i = 0; i < N; ++i) {
    const auto& point = points[i];
    const auto& p = point.basicVector();
    x[i] = p.x();
    y[i] = p.y();
    z[i] = sqr(p.x()) + sqr(p.y());

    // TODO: The following accounts only the hit uncertainty in phi.
    // Albeit it "works nicely", should we account also the
    // uncertainty in r in FPix (and the correlation between phi and r)?
    const float phiErr2 = errors[i].phierr(point);
    const float w = 1.f / (point.perp2() * phiErr2);
    weight[i] = w;

    PRINT << " point " << p.x() << " " << p.y() << " transformed " << x[i] << " " << y[i] << " " << z[i] << " weight "
          << w << std::endl;
  }
  const float invTotWeight = 1.f / std::accumulate(weight.begin(), weight.end(), 0.f);
  PRINT << " invTotWeight " << invTotWeight;

  Eigen::Vector3f mean = Eigen::Vector3f::Zero();
  for (size_t i = 0; i < N; ++i) {
    const float w = weight[i];
    mean[0] += w * x[i];
    mean[1] += w * y[i];
    mean[2] += w * z[i];
  }
  mean *= invTotWeight;
  PRINT << " mean " << mean[0] << " " << mean[1] << " " << mean[2] << std::endl;

  Eigen::Matrix3f A = Eigen::Matrix3f::Zero();
  for (size_t i = 0; i < N; ++i) {
    const auto diff_x = x[i] - mean[0];
    const auto diff_y = y[i] - mean[1];
    const auto diff_z = z[i] - mean[2];
    const auto w = weight[i];

    // exploit the fact that only lower triangular part of the matrix
    // is used in the eigenvalue calculation
    A(0, 0) += w * diff_x * diff_x;
    A(1, 0) += w * diff_y * diff_x;
    A(1, 1) += w * diff_y * diff_y;
    A(2, 0) += w * diff_z * diff_x;
    A(2, 1) += w * diff_z * diff_y;
    A(2, 2) += w * diff_z * diff_z;
    PRINT << " i " << i << " A " << A(0, 0) << " " << A(0, 1) << " " << A(0, 2) << std::endl
          << "     " << A(1, 0) << " " << A(1, 1) << " " << A(1, 2) << std::endl
          << "     " << A(2, 0) << " " << A(2, 1) << " " << A(2, 2) << std::endl;
  }
  A *= 1. / static_cast<float>(N);

  PRINT << " A " << A(0, 0) << " " << A(0, 1) << " " << A(0, 2) << std::endl
        << "   " << A(1, 0) << " " << A(1, 1) << " " << A(1, 2) << std::endl
        << "   " << A(2, 0) << " " << A(2, 1) << " " << A(2, 2) << std::endl;

  // find eigenvalues and vectors
  Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f> eigen(A);
  const auto& eigenValues = eigen.eigenvalues();
  const auto& eigenVectors = eigen.eigenvectors();

  // in Eigen the first eigenvalue is the smallest
  PRINT << " eigenvalues " << eigenValues[0] << " " << eigenValues[1] << " " << eigenValues[2] << std::endl;

  PRINT << " eigen " << eigenVectors(0, 0) << " " << eigenVectors(0, 1) << " " << eigenVectors(0, 2) << std::endl
        << "       " << eigenVectors(1, 0) << " " << eigenVectors(1, 1) << " " << eigenVectors(1, 2) << std::endl
        << "       " << eigenVectors(2, 0) << " " << eigenVectors(2, 1) << " " << eigenVectors(2, 2) << std::endl;

  // eivenvector corresponding smallest eigenvalue
  auto n = eigenVectors.col(0);
  PRINT << " n (unit) " << n[0] << " " << n[1] << " " << n[2] << std::endl;

  const float c = -1.f * (n[0] * mean[0] + n[1] * mean[1] + n[2] * mean[2]);

  PRINT << " c " << c << std::endl;

  // calculate fit parameters
  const auto tmp = 0.5f * 1.f / n[2];
  x0_ = -n[0] * tmp;
  y0_ = -n[1] * tmp;
  const float rho2 = (1 - sqr(n[2]) - 4 * c * n[2]) * sqr(tmp);
  rho_ = std::sqrt(rho2);

  // calculate square sum
  float S = 0;
  for (size_t i = 0; i < N; ++i) {
    const float incr = sqr(c + n[0] * x[i] + n[1] * y[i] + n[2] * z[i]) * weight[i];
#if defined(MK_DEBUG) || defined(EDM_ML_DEBUG)
    const float distance = std::sqrt(sqr(x[i] - x0_) + sqr(y[i] - y0_)) - rho_;
    PRINT << " i " << i << " chi2 incr " << incr << " d(hit, circle) " << distance << " sigma "
          << 1.f / std::sqrt(weight[i]) << std::endl;
#endif
    S += incr;
  }
  chi2_ = S;
}

float FastCircleFit::chi2 ( void  ) const

inline

Definition at line 66 of file FastCircleFit.h.

References chi2_.

Referenced by CAHitQuadrupletGenerator::hitNtuplets().

{ return chi2_; }

float FastCircleFit::rho ( ) const

inline

Definition at line 63 of file FastCircleFit.h.

References rho_.

Referenced by Lepton.Lepton::absIsoFromEA(), and Muon.Muon::absIsoWithFSR().

{ return rho_; }

template<typename T >

T FastCircleFit::sqr ( T  t )

inlineprivate

Definition at line 73 of file FastCircleFit.h.

References submitPVValidationJobs::t.

Referenced by calculate().

             {
    return t * t;
  }

float FastCircleFit::x0 ( ) const

inline

Definition at line 61 of file FastCircleFit.h.

References x0_.

{ return x0_; }

float FastCircleFit::y0 ( ) const

inline

Definition at line 62 of file FastCircleFit.h.

References y0_.

{ return y0_; }

Member Data Documentation

float FastCircleFit::chi2_

private

Definition at line 80 of file FastCircleFit.h.

Referenced by calculate(), and chi2().

float FastCircleFit::rho_

private

Definition at line 79 of file FastCircleFit.h.

Referenced by calculate(), and rho().

float FastCircleFit::x0_

private

Definition at line 77 of file FastCircleFit.h.

Referenced by calculate(), and x0().

float FastCircleFit::y0_

private

Definition at line 78 of file FastCircleFit.h.

Referenced by calculate(), and y0().