File indexing completed on 2024-04-06 12:28:27
0001 #include "RecoTracker/NuclearSeedGenerator/interface/TangentCircle.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #define PI 3.1415926
0005
0006
0007 TangentCircle::TangentCircle(const GlobalVector& direction, const GlobalPoint& inner, const GlobalPoint& outer)
0008 : theInnerPoint(inner), theOuterPoint(outer), theVertexPoint(inner) {
0009 if (theInnerPoint.perp2() > theOuterPoint.perp2()) {
0010 valid = false;
0011 } else
0012 valid = true;
0013
0014 double x1 = inner.x();
0015 double y1 = inner.y();
0016 double x2 = outer.x();
0017 double y2 = outer.y();
0018 double alpha1 = (direction.y() != 0) ? atan(-direction.x() / direction.y()) : PI / 2;
0019 double denominator = 2 * ((x1 - x2) * cos(alpha1) + (y1 - y2) * sin(alpha1));
0020 theRho = (denominator != 0) ? ((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)) / denominator : 1E12;
0021
0022
0023 theX0 = 1E10;
0024 theY0 = 1E10;
0025
0026 theDirectionAtVertex = direction;
0027 theDirectionAtVertex /= theDirectionAtVertex.mag();
0028
0029
0030
0031 theCharge = 0;
0032 theRho = fabs(theRho);
0033
0034 theVertexError = (theInnerPoint - theOuterPoint).mag();
0035 }
0036
0037 TangentCircle::TangentCircle(const GlobalPoint& outerPoint,
0038 const GlobalPoint& innerPoint,
0039 const GlobalPoint& vertexPoint)
0040 : theInnerPoint(innerPoint), theOuterPoint(outerPoint), theVertexPoint(vertexPoint) {
0041 FastCircle circle(outerPoint, innerPoint, vertexPoint);
0042 theX0 = circle.x0();
0043 theY0 = circle.y0();
0044 theRho = circle.rho();
0045 theVertexError = 0;
0046 theCharge = 0;
0047 theDirectionAtVertex = GlobalVector(1000, 1000, 1000);
0048 if (theInnerPoint.perp2() > theOuterPoint.perp2() || !circle.isValid()) {
0049 valid = false;
0050 } else
0051 valid = true;
0052 }
0053
0054 TangentCircle::TangentCircle(const TangentCircle& primCircle,
0055 const GlobalPoint& outerPoint,
0056 const GlobalPoint& innerPoint) {
0057 if (theInnerPoint.perp2() > theOuterPoint.perp2()) {
0058 valid = false;
0059 } else
0060 valid = true;
0061
0062 int NITER = 10;
0063
0064
0065 GlobalPoint InitialVertex(primCircle.outerPoint().x(), primCircle.outerPoint().y(), 0);
0066 GlobalPoint SecInnerPoint(innerPoint.x(), innerPoint.y(), 0);
0067 GlobalPoint SecOuterPoint(outerPoint.x(), outerPoint.y(), 0);
0068
0069
0070 double s = (SecInnerPoint - InitialVertex).mag();
0071 double deltaTheta = s / primCircle.rho();
0072
0073 double minTangentCondition = 1E12;
0074 TangentCircle theCorrectSecCircle;
0075 GlobalPoint vertex = InitialVertex;
0076 int dir = 1;
0077 double theta = deltaTheta / (NITER - 1);
0078
0079 for (int i = 0; i < NITER; i++) {
0080
0081 TangentCircle secCircle(SecOuterPoint, SecInnerPoint, vertex);
0082
0083
0084 double minCond = isTangent(primCircle, secCircle);
0085
0086
0087
0088
0089 if (minCond < minTangentCondition) {
0090 minTangentCondition = minCond;
0091 theCorrectSecCircle = secCircle;
0092 vertex = getPosition(primCircle, secCircle.vertexPoint(), theta, dir);
0093 if (i == 0 && ((vertex - SecInnerPoint).mag() > (InitialVertex - SecInnerPoint).mag())) {
0094 dir = -1;
0095 vertex = getPosition(primCircle, InitialVertex, theta, dir);
0096 LogDebug("NuclearSeedGenerator") << "Change direction to look for vertex"
0097 << "\n";
0098 }
0099 } else
0100 break;
0101 }
0102 theInnerPoint = theCorrectSecCircle.innerPoint();
0103 theOuterPoint = theCorrectSecCircle.outerPoint();
0104 theVertexPoint = theCorrectSecCircle.vertexPoint();
0105 theX0 = theCorrectSecCircle.x0();
0106 theY0 = theCorrectSecCircle.y0();
0107 theRho = theCorrectSecCircle.rho();
0108 theCharge = 0;
0109 theDirectionAtVertex = GlobalVector(1000, 1000, 1000);
0110
0111 theVertexError = s / NITER;
0112 }
0113
0114 double TangentCircle::isTangent(const TangentCircle& primCircle, const TangentCircle& secCircle) const {
0115
0116
0117 double distanceBetweenCircle = (primCircle.x0() - secCircle.x0()) * (primCircle.x0() - secCircle.x0()) +
0118 (primCircle.y0() - secCircle.y0()) * (primCircle.y0() - secCircle.y0());
0119 double RadiusSum = (primCircle.rho() + secCircle.rho()) * (primCircle.rho() + secCircle.rho());
0120 double RadiusDifference = (primCircle.rho() - secCircle.rho()) * (primCircle.rho() - secCircle.rho());
0121
0122 return std::min(fabs(RadiusSum - distanceBetweenCircle), fabs(RadiusDifference - distanceBetweenCircle));
0123 }
0124
0125 GlobalVector TangentCircle::direction(const GlobalPoint& point) const {
0126 if (theY0 > 1E9 || theX0 > 1E9) {
0127 LogDebug("NuclearSeedGenerator") << "Center of TangentCircle not calculated but used !!!"
0128 << "\n";
0129 }
0130
0131
0132 GlobalVector dir(point.y() - theY0, theX0 - point.x(), 0);
0133
0134 dir /= dir.mag();
0135
0136
0137 GlobalVector fastDir = theOuterPoint - theInnerPoint;
0138 double diff = (dir - fastDir).mag();
0139 double sum = (dir + fastDir).mag();
0140
0141 if (sum < diff)
0142 dir = (-1) * dir;
0143
0144 return dir;
0145 }
0146
0147 GlobalVector TangentCircle::directionAtVertex() {
0148 if (theDirectionAtVertex.x() > 999)
0149 theDirectionAtVertex = direction(theVertexPoint);
0150 return theDirectionAtVertex;
0151 }
0152
0153 GlobalPoint TangentCircle::getPosition(const TangentCircle& circle,
0154 const GlobalPoint& initalPosition,
0155 double theta,
0156 int dir) const {
0157 int sign[3];
0158 double x2 = initalPosition.x();
0159 double y2 = initalPosition.y();
0160
0161 if ((x2 > circle.x0()) && dir > 0) {
0162 sign[0] = 1;
0163 sign[1] = -1;
0164 sign[2] = -1;
0165 }
0166 if ((x2 > circle.x0()) && dir < 0) {
0167 sign[0] = 1;
0168 sign[1] = 1;
0169 sign[2] = 1;
0170 }
0171 if ((x2 < circle.x0()) && dir > 0) {
0172 sign[0] = -1;
0173 sign[1] = 1;
0174 sign[2] = -1;
0175 }
0176 if ((x2 < circle.x0()) && dir < 0) {
0177 sign[0] = -1;
0178 sign[1] = -1;
0179 sign[2] = 1;
0180 }
0181
0182 double l = 2 * circle.rho() * sin(theta / 2);
0183 double alpha = atan((y2 - circle.y0()) / (x2 - circle.x0()));
0184 double beta = PI / 2 - theta / 2;
0185 double gamma = PI + sign[2] * alpha - beta;
0186
0187 double xnew = x2 + sign[0] * l * cos(gamma);
0188 double ynew = y2 + sign[1] * l * sin(gamma);
0189
0190 return GlobalPoint(xnew, ynew, 0);
0191 }
0192
0193 double TangentCircle::curvatureError() {
0194 if ((theInnerPoint - theVertexPoint).mag() < theVertexError) {
0195 TangentCircle circle1(directionAtVertex(), theVertexPoint - theVertexError * directionAtVertex(), theOuterPoint);
0196 TangentCircle circle2(directionAtVertex(), theVertexPoint + theVertexError * directionAtVertex(), theOuterPoint);
0197 return fabs(1 / circle1.rho() - 1 / circle2.rho());
0198 } else {
0199 TangentCircle circle1(theOuterPoint, theInnerPoint, theVertexPoint - theVertexError * directionAtVertex());
0200 TangentCircle circle2(theOuterPoint, theInnerPoint, theVertexPoint + theVertexError * directionAtVertex());
0201 return fabs(1 / circle1.rho() - 1 / circle2.rho());
0202 }
0203 }
0204
0205 int TangentCircle::charge(float magz) {
0206 if (theCharge != 0)
0207 return theCharge;
0208
0209 if (theX0 > 1E9 || theY0 > 1E9)
0210 theCharge = chargeLocally(magz, directionAtVertex());
0211 else {
0212 GlobalPoint center(theX0, theY0, 0);
0213 GlobalVector u = center - theVertexPoint;
0214 GlobalVector v = directionAtVertex();
0215
0216
0217 GlobalVector F(v.y() * magz, -v.x() * magz, 0);
0218 if (u.x() * F.x() + u.y() * F.y() > 0)
0219 theCharge = -1;
0220 else
0221 theCharge = 1;
0222
0223 if (theCharge != chargeLocally(magz, v)) {
0224 LogDebug("NuclearSeedGenerator") << "Inconsistency in calculation of the charge"
0225 << "\n";
0226 }
0227 }
0228 return theCharge;
0229 }
0230
0231 int TangentCircle::chargeLocally(float magz, GlobalVector v) const {
0232 GlobalVector u = theOuterPoint - theVertexPoint;
0233 double tz = v.x() * u.y() - v.y() * u.x();
0234
0235 if (tz * magz > 0)
0236 return 1;
0237 else
0238 return -1;
0239 }