PixelSeeding/plugins/ThirdHitPredictionFromInvLine.cc

0001 #include "ThirdHitPredictionFromInvLine.h"
0002
0003 #include <cmath>
0004 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0005 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0006
0007 #include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
0008 #include "RecoTracker/TkMSParametrization/interface/PixelRecoRange.h"
0009
0010 template <class T>
0011 T sqr(T t) {
0012   return t * t;
0013 }
0014
0015 typedef Basic3DVector<double> Point3D;
0016 typedef Basic2DVector<double> Point2D;
0017 typedef PixelRecoRange<double> Ranged;
0018
0019 #include <iostream>
0020
0021 using namespace std;
0022
0023 ThirdHitPredictionFromInvLine::ThirdHitPredictionFromInvLine(const GlobalPoint &P1,
0024                                                              const GlobalPoint &P2,
0025                                                              double errorRPhiP1,
0026                                                              double errorRPhiP2)
0027     : nPoints(0),
0028       theSum(0.),
0029       theSumU(0.),
0030       theSumUU(0.),
0031       theSumV(0.),
0032       theSumUV(0.),
0033       theSumVV(0.),
0034       hasParameters(false),
0035       theCurvatureValue(0.),
0036       theCurvatureError(0.),
0037       theChi2(0.) {
0038   GlobalVector aX = GlobalVector(P1.x(), P1.y(), 0.).unit();
0039   GlobalVector aY(-aX.y(), aX.x(), 0.);
0040   GlobalVector aZ(0., 0., 1.);
0041   theRotation = Rotation(aX, aY, aZ);
0042
0043   add(P1, errorRPhiP1);
0044   add(P2, errorRPhiP2);
0045 }
0046
0047 GlobalPoint ThirdHitPredictionFromInvLine::crossing(double radius) const {
0048   double A = -(theSum * theSumUV - theSumU * theSumV) / (sqr(theSumU) - theSum * theSumUU);
0049   double B = (theSumU * theSumUV - theSumUU * theSumV) / (sqr(theSumU) - theSum * theSumUU);
0050   double delta = sqr(2. * A * B) - 4 * (1 + sqr(A)) * (sqr(B) - sqr(1 / radius));
0051   double sqrtdelta = (delta > 0.) ? sqrt(delta) : 0.;
0052   double u1 = (-2. * A * B + sqrtdelta) / 2. / (1 + sqr(A));
0053   double v1 = A * u1 + B;
0054   Point2D tmp = PointUV(u1, v1, &theRotation).unmap();
0055   return GlobalPoint(tmp.x(), tmp.y(), 0.);
0056 }
0057
0058 void ThirdHitPredictionFromInvLine::add(const GlobalPoint &p, double errorRPhi) {
0059   double weigth = sqr(sqr(p.perp()) / errorRPhi);
0060   add(PointUV(Point2D(p.x(), p.y()), &theRotation), weigth);
0061 }
0062
0063 void ThirdHitPredictionFromInvLine::add(const ThirdHitPredictionFromInvLine::PointUV &point, double weigth) {
0064   hasParameters = false;
0065   nPoints++;
0066   theSum += weigth;
0067   theSumU += point.u() * weigth;
0068   theSumUU += sqr(point.u()) * weigth;
0069   theSumV += point.v() * weigth;
0070   theSumUV += point.u() * point.v() * weigth;
0071   theSumVV += sqr(point.v()) * weigth;
0072   check();
0073 }
0074
0075 void ThirdHitPredictionFromInvLine::remove(const GlobalPoint &p, double errorRPhi) {
0076   hasParameters = false;
0077   PointUV point(Point2D(p.x(), p.y()), &theRotation);
0078   double weigth = sqr(sqr(p.perp()) / errorRPhi);
0079   nPoints--;
0080   theSum -= weigth;
0081   theSumU -= point.u() * weigth;
0082   theSumUU -= sqr(point.u()) * weigth;
0083   theSumV -= point.v() * weigth;
0084   theSumUV -= point.u() * point.v() * weigth;
0085   theSumVV -= sqr(point.v()) * weigth;
0086   check();
0087 }
0088
0089 void ThirdHitPredictionFromInvLine::print() const {
0090   std::cout << " nPoints: " << nPoints << " theSumU: " << theSumU << " theSumUU: " << theSumUU
0091             << " theSumV: " << theSumV << " theSumUV: " << theSumUV << std::endl;
0092 }
0093
0094 void ThirdHitPredictionFromInvLine::check() {
0095   if (hasParameters)
0096     return;
0097
0098   long double D = theSumUU * theSum - theSumU * theSumU;
0099   long double A = (theSumUV * theSum - theSumU * theSumV) / D;
0100   long double B = (theSumUU * theSumV - theSumUV * theSumU) / D;
0101   double rho = 2. * fabs(B) / sqrt(1 + sqr(A));
0102   double sigmaA2 = theSum / D;
0103   double sigmaB2 = theSumUU / D;
0104
0105   hasParameters = true;
0106   theCurvatureError = sqrt(sqr(rho / B) * sigmaB2 + sqr(rho / (1 + sqr(A))) * sigmaA2);
0107   theCurvatureValue = 2. * fabs(B) / sqrt(1 + sqr(A));
0108   theChi2 = theSumVV - 2 * A * theSumUV - 2 * B * theSumV + 2 * A * B * theSumU + B * B * theSum + A * A * theSumUU;
0109 }
0110
0111 /*
0112 GlobalPoint ThirdHitPredictionFromInvLine::center() const
0113 {
0114   long double den = theSumU*theSumUV - theSumUU*theSumV;
0115   double a = (theSum*theSumUV-theSumU*theSumV)/2./den;
0116   double b = (sqr(theSumU)-theSum*theSumUU)/2./den;
0117   Point3D tmp = theRotation.multiplyInverse( Point2D(a,b) );
0118   return GlobalPoint(tmp.x(), tmp.y(), 0.);
0119 }
0120 */