ThirdHitPredictionFromInvLine.cc

Go to the documentation of this file.

#include "ThirdHitPredictionFromInvLine.h"

#include <cmath>
#include "DataFormats/GeometryVector/interface/GlobalVector.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"

#include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
#include "RecoTracker/TkMSParametrization/interface/PixelRecoRange.h"

template <class T>
T sqr(T t) {
  return t * t;
}

typedef Basic3DVector<double> Point3D;
typedef Basic2DVector<double> Point2D;
typedef PixelRecoRange<double> Ranged;

#include <iostream>

using namespace std;

ThirdHitPredictionFromInvLine::ThirdHitPredictionFromInvLine(const GlobalPoint &P1,
                                                             const GlobalPoint &P2,
                                                             double errorRPhiP1,
                                                             double errorRPhiP2)
    : nPoints(0),
      theSum(0.),
      theSumU(0.),
      theSumUU(0.),
      theSumV(0.),
      theSumUV(0.),
      theSumVV(0.),
      hasParameters(false),
      theCurvatureValue(0.),
      theCurvatureError(0.),
      theChi2(0.) {
  GlobalVector aX = GlobalVector(P1.x(), P1.y(), 0.).unit();
  GlobalVector aY(-aX.y(), aX.x(), 0.);
  GlobalVector aZ(0., 0., 1.);
  theRotation = Rotation(aX, aY, aZ);

  add(P1, errorRPhiP1);
  add(P2, errorRPhiP2);
}

GlobalPoint ThirdHitPredictionFromInvLine::crossing(double radius) const {
  double A = -(theSum * theSumUV - theSumU * theSumV) / (sqr(theSumU) - theSum * theSumUU);
  double B = (theSumU * theSumUV - theSumUU * theSumV) / (sqr(theSumU) - theSum * theSumUU);
  double delta = sqr(2. * A * B) - 4 * (1 + sqr(A)) * (sqr(B) - sqr(1 / radius));
  double sqrtdelta = (delta > 0.) ? sqrt(delta) : 0.;
  double u1 = (-2. * A * B + sqrtdelta) / 2. / (1 + sqr(A));
  double v1 = A * u1 + B;
  Point2D tmp = PointUV(u1, v1, &theRotation).unmap();
  return GlobalPoint(tmp.x(), tmp.y(), 0.);
}

void ThirdHitPredictionFromInvLine::add(const GlobalPoint &p, double errorRPhi) {
  double weigth = sqr(sqr(p.perp()) / errorRPhi);
  add(PointUV(Point2D(p.x(), p.y()), &theRotation), weigth);
}

void ThirdHitPredictionFromInvLine::add(const ThirdHitPredictionFromInvLine::PointUV &point, double weigth) {
  hasParameters = false;
  nPoints++;
  theSum += weigth;
  theSumU += point.u() * weigth;
  theSumUU += sqr(point.u()) * weigth;
  theSumV += point.v() * weigth;
  theSumUV += point.u() * point.v() * weigth;
  theSumVV += sqr(point.v()) * weigth;
  check();
}

void ThirdHitPredictionFromInvLine::remove(const GlobalPoint &p, double errorRPhi) {
  hasParameters = false;
  PointUV point(Point2D(p.x(), p.y()), &theRotation);
  double weigth = sqr(sqr(p.perp()) / errorRPhi);
  nPoints--;
  theSum -= weigth;
  theSumU -= point.u() * weigth;
  theSumUU -= sqr(point.u()) * weigth;
  theSumV -= point.v() * weigth;
  theSumUV -= point.u() * point.v() * weigth;
  theSumVV -= sqr(point.v()) * weigth;
  check();
}

void ThirdHitPredictionFromInvLine::print() const {
  std::cout << " nPoints: " << nPoints << " theSumU: " << theSumU << " theSumUU: " << theSumUU
            << " theSumV: " << theSumV << " theSumUV: " << theSumUV << std::endl;
}

void ThirdHitPredictionFromInvLine::check() {
  if (hasParameters)
    return;

  long double D = theSumUU * theSum - theSumU * theSumU;
  long double A = (theSumUV * theSum - theSumU * theSumV) / D;
  long double B = (theSumUU * theSumV - theSumUV * theSumU) / D;
  double rho = 2. * fabs(B) / sqrt(1 + sqr(A));
  double sigmaA2 = theSum / D;
  double sigmaB2 = theSumUU / D;

  hasParameters = true;
  theCurvatureError = sqrt(sqr(rho / B) * sigmaB2 + sqr(rho / (1 + sqr(A))) * sigmaA2);
  theCurvatureValue = 2. * fabs(B) / sqrt(1 + sqr(A));
  theChi2 = theSumVV - 2 * A * theSumUV - 2 * B * theSumV + 2 * A * B * theSumU + B * B * theSum + A * A * theSumUU;
}

/*
GlobalPoint ThirdHitPredictionFromInvLine::center() const
{
  long double den = theSumU*theSumUV - theSumUU*theSumV; 
  double a = (theSum*theSumUV-theSumU*theSumV)/2./den;
  double b = (sqr(theSumU)-theSum*theSumUU)/2./den;
  Point3D tmp = theRotation.multiplyInverse( Point2D(a,b) );
  return GlobalPoint(tmp.x(), tmp.y(), 0.);
}
*/