NuclearSeedGenerator/src/TangentCircle.cc

0001 #include "RecoTracker/NuclearSeedGenerator/interface/TangentCircle.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #define PI 3.1415926
0005
0006 // TODO: is not valid don't do any calculations and return init values
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   // TODO : variable not yet calculated look in nucl.C
0023   theX0 = 1E10;
0024   theY0 = 1E10;
0025
0026   theDirectionAtVertex = direction;
0027   theDirectionAtVertex /= theDirectionAtVertex.mag();
0028
0029   //theCharge = (theRho>0) ? -1 : 1;
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   // Initial vertex used = outerPoint of the primary circle (should be the first estimation of the nuclear interaction position)
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   // distance between the initial vertex and the inner point of the secondary circle
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     // get the circle which pass through outerPoint, innerPoint and the vertex
0081     TangentCircle secCircle(SecOuterPoint, SecInnerPoint, vertex);
0082
0083     // get a value relative to the tangentness of the 2 circles
0084     double minCond = isTangent(primCircle, secCircle);
0085
0086     // double dirDiff = (primCircle.direction(vertex) - secCircle.direction(vertex)).mag();
0087     // if( dirDiff > 1) dirDiff = 2-dirDiff;
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   // return a value that should be equal to 0 if primCircle and secCircle are tangent
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   // calculate the direction perpendicular to the vector v = point - center_of_circle
0132   GlobalVector dir(point.y() - theY0, theX0 - point.x(), 0);
0133
0134   dir /= dir.mag();
0135
0136   // Check the sign :
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     // F = force vector
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 }