TangentCircle.cc

Go to the documentation of this file.

#include "RecoTracker/NuclearSeedGenerator/interface/TangentCircle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#define PI 3.1415926

// TODO: is not valid don't do any calculations and return init values
TangentCircle::TangentCircle(const GlobalVector& direction, const GlobalPoint& inner, const GlobalPoint& outer)
    : theInnerPoint(inner), theOuterPoint(outer), theVertexPoint(inner) {
  if (theInnerPoint.perp2() > theOuterPoint.perp2()) {
    valid = false;
  } else
    valid = true;

  double x1 = inner.x();
  double y1 = inner.y();
  double x2 = outer.x();
  double y2 = outer.y();
  double alpha1 = (direction.y() != 0) ? atan(-direction.x() / direction.y()) : PI / 2;
  double denominator = 2 * ((x1 - x2) * cos(alpha1) + (y1 - y2) * sin(alpha1));
  theRho = (denominator != 0) ? ((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)) / denominator : 1E12;

  // TODO : variable not yet calculated look in nucl.C
  theX0 = 1E10;
  theY0 = 1E10;

  theDirectionAtVertex = direction;
  theDirectionAtVertex /= theDirectionAtVertex.mag();

  //theCharge = (theRho>0) ? -1 : 1;

  theCharge = 0;
  theRho = fabs(theRho);

  theVertexError = (theInnerPoint - theOuterPoint).mag();
}

TangentCircle::TangentCircle(const GlobalPoint& outerPoint,
                             const GlobalPoint& innerPoint,
                             const GlobalPoint& vertexPoint)
    : theInnerPoint(innerPoint), theOuterPoint(outerPoint), theVertexPoint(vertexPoint) {
  FastCircle circle(outerPoint, innerPoint, vertexPoint);
  theX0 = circle.x0();
  theY0 = circle.y0();
  theRho = circle.rho();
  theVertexError = 0;
  theCharge = 0;
  theDirectionAtVertex = GlobalVector(1000, 1000, 1000);
  if (theInnerPoint.perp2() > theOuterPoint.perp2() || !circle.isValid()) {
    valid = false;
  } else
    valid = true;
}

TangentCircle::TangentCircle(const TangentCircle& primCircle,
                             const GlobalPoint& outerPoint,
                             const GlobalPoint& innerPoint) {
  if (theInnerPoint.perp2() > theOuterPoint.perp2()) {
    valid = false;
  } else
    valid = true;

  int NITER = 10;

  // Initial vertex used = outerPoint of the primary circle (should be the first estimation of the nuclear interaction position)
  GlobalPoint InitialVertex(primCircle.outerPoint().x(), primCircle.outerPoint().y(), 0);
  GlobalPoint SecInnerPoint(innerPoint.x(), innerPoint.y(), 0);
  GlobalPoint SecOuterPoint(outerPoint.x(), outerPoint.y(), 0);

  // distance between the initial vertex and the inner point of the secondary circle
  double s = (SecInnerPoint - InitialVertex).mag();
  double deltaTheta = s / primCircle.rho();

  double minTangentCondition = 1E12;
  TangentCircle theCorrectSecCircle;
  GlobalPoint vertex = InitialVertex;
  int dir = 1;
  double theta = deltaTheta / (NITER - 1);

  for (int i = 0; i < NITER; i++) {
    // get the circle which pass through outerPoint, innerPoint and the vertex
    TangentCircle secCircle(SecOuterPoint, SecInnerPoint, vertex);

    // get a value relative to the tangentness of the 2 circles
    double minCond = isTangent(primCircle, secCircle);

    // double dirDiff = (primCircle.direction(vertex) - secCircle.direction(vertex)).mag();
    // if( dirDiff > 1) dirDiff = 2-dirDiff;

    if (minCond < minTangentCondition) {
      minTangentCondition = minCond;
      theCorrectSecCircle = secCircle;
      vertex = getPosition(primCircle, secCircle.vertexPoint(), theta, dir);
      if (i == 0 && ((vertex - SecInnerPoint).mag() > (InitialVertex - SecInnerPoint).mag())) {
        dir = -1;
        vertex = getPosition(primCircle, InitialVertex, theta, dir);
        LogDebug("NuclearSeedGenerator") << "Change direction to look for vertex"
                                         << "\n";
      }
    } else
      break;
  }
  theInnerPoint = theCorrectSecCircle.innerPoint();
  theOuterPoint = theCorrectSecCircle.outerPoint();
  theVertexPoint = theCorrectSecCircle.vertexPoint();
  theX0 = theCorrectSecCircle.x0();
  theY0 = theCorrectSecCircle.y0();
  theRho = theCorrectSecCircle.rho();
  theCharge = 0;
  theDirectionAtVertex = GlobalVector(1000, 1000, 1000);

  theVertexError = s / NITER;
}

double TangentCircle::isTangent(const TangentCircle& primCircle, const TangentCircle& secCircle) const {
  // return a value that should be equal to 0 if primCircle and secCircle are tangent

  double distanceBetweenCircle = (primCircle.x0() - secCircle.x0()) * (primCircle.x0() - secCircle.x0()) +
                                 (primCircle.y0() - secCircle.y0()) * (primCircle.y0() - secCircle.y0());
  double RadiusSum = (primCircle.rho() + secCircle.rho()) * (primCircle.rho() + secCircle.rho());
  double RadiusDifference = (primCircle.rho() - secCircle.rho()) * (primCircle.rho() - secCircle.rho());

  return std::min(fabs(RadiusSum - distanceBetweenCircle), fabs(RadiusDifference - distanceBetweenCircle));
}

GlobalVector TangentCircle::direction(const GlobalPoint& point) const {
  if (theY0 > 1E9 || theX0 > 1E9) {
    LogDebug("NuclearSeedGenerator") << "Center of TangentCircle not calculated but used !!!"
                                     << "\n";
  }

  // calculate the direction perpendicular to the vector v = point - center_of_circle
  GlobalVector dir(point.y() - theY0, theX0 - point.x(), 0);

  dir /= dir.mag();

  // Check the sign :
  GlobalVector fastDir = theOuterPoint - theInnerPoint;
  double diff = (dir - fastDir).mag();
  double sum = (dir + fastDir).mag();

  if (sum < diff)
    dir = (-1) * dir;

  return dir;
}

GlobalVector TangentCircle::directionAtVertex() {
  if (theDirectionAtVertex.x() > 999)
    theDirectionAtVertex = direction(theVertexPoint);
  return theDirectionAtVertex;
}

GlobalPoint TangentCircle::getPosition(const TangentCircle& circle,
                                       const GlobalPoint& initalPosition,
                                       double theta,
                                       int dir) const {
  int sign[3];
  double x2 = initalPosition.x();
  double y2 = initalPosition.y();

  if ((x2 > circle.x0()) && dir > 0) {
    sign[0] = 1;
    sign[1] = -1;
    sign[2] = -1;
  }
  if ((x2 > circle.x0()) && dir < 0) {
    sign[0] = 1;
    sign[1] = 1;
    sign[2] = 1;
  }
  if ((x2 < circle.x0()) && dir > 0) {
    sign[0] = -1;
    sign[1] = 1;
    sign[2] = -1;
  }
  if ((x2 < circle.x0()) && dir < 0) {
    sign[0] = -1;
    sign[1] = -1;
    sign[2] = 1;
  }

  double l = 2 * circle.rho() * sin(theta / 2);
  double alpha = atan((y2 - circle.y0()) / (x2 - circle.x0()));
  double beta = PI / 2 - theta / 2;
  double gamma = PI + sign[2] * alpha - beta;

  double xnew = x2 + sign[0] * l * cos(gamma);
  double ynew = y2 + sign[1] * l * sin(gamma);

  return GlobalPoint(xnew, ynew, 0);
}

double TangentCircle::curvatureError() {
  if ((theInnerPoint - theVertexPoint).mag() < theVertexError) {
    TangentCircle circle1(directionAtVertex(), theVertexPoint - theVertexError * directionAtVertex(), theOuterPoint);
    TangentCircle circle2(directionAtVertex(), theVertexPoint + theVertexError * directionAtVertex(), theOuterPoint);
    return fabs(1 / circle1.rho() - 1 / circle2.rho());
  } else {
    TangentCircle circle1(theOuterPoint, theInnerPoint, theVertexPoint - theVertexError * directionAtVertex());
    TangentCircle circle2(theOuterPoint, theInnerPoint, theVertexPoint + theVertexError * directionAtVertex());
    return fabs(1 / circle1.rho() - 1 / circle2.rho());
  }
}

int TangentCircle::charge(float magz) {
  if (theCharge != 0)
    return theCharge;

  if (theX0 > 1E9 || theY0 > 1E9)
    theCharge = chargeLocally(magz, directionAtVertex());
  else {
    GlobalPoint center(theX0, theY0, 0);
    GlobalVector u = center - theVertexPoint;
    GlobalVector v = directionAtVertex();

    // F = force vector
    GlobalVector F(v.y() * magz, -v.x() * magz, 0);
    if (u.x() * F.x() + u.y() * F.y() > 0)
      theCharge = -1;
    else
      theCharge = 1;

    if (theCharge != chargeLocally(magz, v)) {
      LogDebug("NuclearSeedGenerator") << "Inconsistency in calculation of the charge"
                                       << "\n";
    }
  }
  return theCharge;
}

int TangentCircle::chargeLocally(float magz, GlobalVector v) const {
  GlobalVector u = theOuterPoint - theVertexPoint;
  double tz = v.x() * u.y() - v.y() * u.x();

  if (tz * magz > 0)
    return 1;
  else
    return -1;
}