CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/RecoVertex/GhostTrackFitter/src/TrackGhostTrackState.cc

Go to the documentation of this file.
00001 #include <cmath>
00002 
00003 #include <Math/SMatrix.h>
00004 #include <Math/MatrixFunctions.h>
00005 
00006 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00007 #include "DataFormats/GeometryVector/interface/GlobalVector.h" 
00008 #include "DataFormats/GeometrySurface/interface/Line.h"
00009 
00010 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00011 #include "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h"
00012 #include "TrackingTools/GeomPropagators/interface/AnalyticalTrajectoryExtrapolatorToLine.h"
00013 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00014 #include "TrackingTools/TrajectoryParametrization/interface/CartesianTrajectoryError.h"
00015 
00016 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackPrediction.h"
00017 
00018 #include "RecoVertex/GhostTrackFitter/interface/TrackGhostTrackState.h"
00019 
00020 using namespace reco;
00021 
00022 namespace {
00023         static inline double sqr(double arg) { return arg * arg; }
00024 
00025         using namespace ROOT::Math;
00026 
00027         typedef SVector<double, 3> Vector3;
00028 
00029         typedef SMatrix<double, 3, 3, MatRepSym<double, 3> > Matrix3S;
00030         typedef SMatrix<double, 6, 6, MatRepSym<double, 6> > Matrix6S;
00031         typedef SMatrix<double, 3, 3> Matrix33;
00032         typedef SMatrix<double, 3, 6> Matrix36;
00033 
00034         static inline Vector3 conv(const GlobalVector &vec)
00035         {
00036                 Vector3 result;
00037                 result[0] = vec.x();
00038                 result[1] = vec.y();
00039                 result[2] = vec.z();
00040                 return result;
00041         }
00042 }
00043 
00044 bool TrackGhostTrackState::linearize(const GhostTrackPrediction &pred,
00045                                      bool initial, double lambda)
00046 {
00047         AnalyticalTrajectoryExtrapolatorToLine extrap(track_.field());
00048 
00049         GlobalPoint origin = pred.origin();
00050         GlobalVector direction = pred.direction();
00051 
00052         if (isValid() && !initial) {
00053                 GlobalPoint pca = origin + lambda_ * direction;
00054                 Line line(pca, direction);
00055                 tsos_ = extrap.extrapolate(tsos_, line);
00056         } else {
00057                 GlobalPoint pca = origin + lambda * direction;
00058                 Line line(pca, direction);
00059                 tsos_ = extrap.extrapolate(track_.impactPointState(), line);
00060         }
00061 
00062         if (!isValid())
00063                 return false;
00064 
00065         lambda_ = (tsos_.globalPosition() - origin) * direction / pred.rho2();
00066 
00067         return true;
00068 }
00069 
00070 bool TrackGhostTrackState::linearize(const GhostTrackPrediction &pred,
00071                                      double lambda)
00072 {
00073         AnalyticalImpactPointExtrapolator extrap(track_.field());
00074 
00075         GlobalPoint point = pred.position(lambda);
00076 
00077         tsos_ = extrap.extrapolate(track_.impactPointState(), point);
00078         if (!isValid())
00079                 return false;
00080 
00081         lambda_ = lambda;
00082 
00083         return true;
00084 }
00085 
00086 BasicGhostTrackState::Vertex TrackGhostTrackState::vertexStateOnGhostTrack(
00087         const GhostTrackPrediction &pred, bool withMeasurementError) const
00088 {
00089         using namespace ROOT::Math;
00090 
00091         if (!isValid())
00092                 return Vertex();
00093 
00094         GlobalPoint origin = pred.origin();
00095         GlobalVector direction = pred.direction();
00096 
00097         double rho2 = pred.rho2();
00098         double rho = std::sqrt(rho2);
00099         double lambda = (tsos_.globalPosition() - origin) * direction / rho2;
00100         GlobalPoint pos = origin + lambda * direction;
00101 
00102         GlobalVector momentum = tsos_.globalMomentum();
00103         double mom = momentum.mag();
00104 
00105         Vector3 b = conv(direction) / rho;
00106         Vector3 d = conv(momentum) / mom;
00107         double l = Dot(b, d);
00108         double g = 1. / (1. - sqr(l));
00109 
00110         Vector3 ca = conv(pos - tsos_.globalPosition());
00111         Vector3 bd = b - l * d;
00112         b *= g;
00113 
00114         Matrix33 pA = TensorProd(b, bd);
00115         Matrix33 pB = TensorProd(b, ca);
00116 
00117         Matrix36 jacobian;
00118         jacobian.Place_at(-pA + Matrix33(SMatrixIdentity()), 0, 0);
00119         jacobian.Place_at(pB / rho, 0, 3);
00120         Matrix3S error = Similarity(jacobian, pred.cartesianError(lambda));
00121 
00122         if (withMeasurementError) {
00123                 jacobian.Place_at(pA, 0, 0);
00124                 jacobian.Place_at(-pB / mom, 0, 3);
00125                 error += Similarity(jacobian, tsos_.cartesianError().matrix());
00126         }
00127 
00128         return Vertex(pos, error);
00129 }
00130 
00131 BasicGhostTrackState::Vertex TrackGhostTrackState::vertexStateOnMeasurement(
00132         const GhostTrackPrediction &pred, bool withGhostTrackError) const
00133 {
00134         using namespace ROOT::Math;
00135 
00136         if (!isValid())
00137                 return Vertex();
00138 
00139         GlobalPoint origin = pred.origin();
00140         GlobalVector direction = pred.direction();
00141 
00142         double rho2 = pred.rho2();
00143         double rho = std::sqrt(rho2);
00144         double lambda = (tsos_.globalPosition() - origin) * direction / rho2;
00145         GlobalPoint pos = origin + lambda * direction;
00146 
00147         GlobalVector momentum = tsos_.globalMomentum();
00148         double mom = momentum.mag();
00149 
00150         Vector3 b = conv(direction) / rho;
00151         Vector3 d = conv(momentum) / mom;
00152         double l = Dot(b, d);
00153         double g = 1. / (1. - sqr(l));
00154 
00155         Vector3 ca = conv(tsos_.globalPosition() - pos);
00156         Vector3 bd = l * b - d;
00157         d *= g;
00158 
00159         Matrix33 pC = TensorProd(d, bd);
00160         Matrix33 pD = TensorProd(d, ca);
00161 
00162         Matrix36 jacobian;
00163         jacobian.Place_at(pC + Matrix33(SMatrixIdentity()), 0, 0);
00164         jacobian.Place_at(pD / mom, 0, 3);
00165         Matrix3S error = Similarity(jacobian, tsos_.cartesianError().matrix());
00166 
00167         if (withGhostTrackError) {
00168                 jacobian.Place_at(-pC, 0, 0);
00169                 jacobian.Place_at(-pD / rho, 0, 3);
00170                 error += Similarity(jacobian, pred.cartesianError(lambda));
00171         }
00172 
00173         return Vertex(tsos_.globalPosition(), error);
00174 }