BTagTools/src/SignedImpactParameter3D.cc

0001 #include "DataFormats/GeometrySurface/interface/Line.h"
0002
0003 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0004 #include "TrackingTools/GeomPropagators/interface/AnalyticalTrajectoryExtrapolatorToLine.h"
0005 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0006 #include "RecoBTag/BTagTools/interface/SignedImpactParameter3D.h"
0007 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0008
0009 #include "CLHEP/Vector/ThreeVector.h"
0010 #include "CLHEP/Vector/LorentzVector.h"
0011 #include "CLHEP/Matrix/Vector.h"
0012 #include <string>
0013
0014 using namespace std;
0015 using namespace reco;
0016
0017 pair<bool, Measurement1D> SignedImpactParameter3D::apply(const TransientTrack& transientTrack,
0018                                                          const GlobalVector& direction,
0019                                                          const Vertex& vertex) const {
0020   double theValue = 0.;
0021   double theError = 0.;
0022   bool theIsValid = false;
0023
0024   TrajectoryStateOnSurface TSOS = transientTrack.impactPointState();
0025
0026   if (!TSOS.isValid()) {
0027     cout << "====>>>> SignedImpactParameter3D::apply : TSOS not valid = " << TSOS.isValid() << endl;
0028     return pair<bool, Measurement1D>(theIsValid, Measurement1D(0., 0.));
0029   }
0030
0031   const FreeTrajectoryState* FTS = TSOS.freeTrajectoryState();
0032
0033   const GlobalVector& JetDirection(direction);
0034
0035   TrajectoryStateOnSurface theTSOS = closestApproachToJet(*FTS, vertex, JetDirection, transientTrack.field());
0036   theIsValid = theTSOS.isValid();
0037
0038   if (theIsValid) {
0039     GlobalVector D = distance(theTSOS, vertex, JetDirection);
0040     GlobalVector J = JetDirection.unit();
0041     GlobalPoint vertexPosition(vertex.x(), vertex.y(), vertex.z());
0042     double theDistanceAlongJetAxis = J.dot(theTSOS.globalPosition() - vertexPosition);
0043     theValue = D.mag() * (theDistanceAlongJetAxis / abs(theDistanceAlongJetAxis));
0044
0045     GlobalVector DD = D.unit();
0046     GlobalPoint T0 = theTSOS.globalPosition();
0047     GlobalVector T1 = theTSOS.globalMomentum();
0048     GlobalVector TT1 = T1.unit();
0049     GlobalVector Xi(T0.x() - vertex.position().x(), T0.y() - vertex.position().y(), T0.z() - vertex.position().z());
0050
0051     AlgebraicVector6 deriv;
0052     AlgebraicVector3 deriv_v;
0053
0054     deriv_v[0] = -DD.x();
0055     deriv_v[1] = -DD.y();
0056     deriv_v[2] = -DD.z();
0057
0058     deriv[0] = DD.x();
0059     deriv[1] = DD.y();
0060     deriv[2] = DD.z();
0061     deriv[3] = -(TT1.dot(Xi) * DD.x()) / T1.mag();
0062     deriv[4] = -(TT1.dot(Xi) * DD.y()) / T1.mag();
0063     deriv[5] = -(TT1.dot(Xi) * DD.z()) / T1.mag();
0064
0065     double E1 = ROOT::Math::Similarity(deriv, theTSOS.cartesianError().matrix());
0066     double E2 = ROOT::Math::Similarity(deriv_v, vertex.covariance());
0067     //    double E2 = RecoVertex::convertError(vertex.covariance()).matrix().similarity(deriv_v);
0068     //    double E2 = 0.; // no vertex error because of stupid use of hundreds of different types for same thing
0069     theError = sqrt(E1 + E2);
0070
0071     Measurement1D A(theValue, theError);
0072
0073     return pair<bool, Measurement1D>(theIsValid, A);
0074   } else {
0075     return pair<bool, Measurement1D>(theIsValid, Measurement1D(0., 0.));
0076   }  // endif (isValid)
0077 }
0078
0079 TrajectoryStateOnSurface SignedImpactParameter3D::closestApproachToJet(const FreeTrajectoryState& aFTS,
0080                                                                        const Vertex& vertex,
0081                                                                        const GlobalVector& aJetDirection,
0082                                                                        const MagneticField* field) {
0083   GlobalVector J = aJetDirection.unit();
0084
0085   Line::PositionType pos(GlobalPoint(vertex.x(), vertex.y(), vertex.z()));
0086   Line::DirectionType dir(J);
0087   Line Jet(pos, dir);
0088
0089   AnalyticalTrajectoryExtrapolatorToLine TETL(field);
0090
0091   return TETL.extrapolate(aFTS, Jet);
0092 }
0093
0094 GlobalVector SignedImpactParameter3D::distance(const TrajectoryStateOnSurface& aTSOS,
0095                                                const Vertex& vertex,
0096                                                const GlobalVector& aJetDirection) {
0097   Line::PositionType pos2(aTSOS.globalPosition());
0098   Line::DirectionType dir2((aTSOS.globalMomentum()).unit());
0099   Line T(pos2, dir2);
0100
0101   GlobalPoint X = GlobalPoint(vertex.x(), vertex.y(), vertex.z());  // aVertex.position();
0102
0103   GlobalVector D = T.distance(X);
0104
0105   return D;
0106 }
0107
0108 pair<double, Measurement1D> SignedImpactParameter3D::distanceWithJetAxis(const TransientTrack& track,
0109                                                                          const GlobalVector& direction,
0110                                                                          const Vertex& vertex) {
0111   double theDistanceAlongJetAxis(0.);
0112   double theDistanceToJetAxis(0.);
0113   double theLDist_err(0.);
0114   TrajectoryStateOnSurface TSOS = track.impactPointState();
0115
0116   if (!TSOS.isValid()) {
0117     cout << "====>>>> SignedImpactParameter3D::distanceWithJetAxis : TSOS not valid = " << TSOS.isValid() << endl;
0118     return pair<double, Measurement1D>(theDistanceAlongJetAxis, Measurement1D(theDistanceToJetAxis, theLDist_err));
0119   }
0120
0121   const FreeTrajectoryState* FTS = TSOS.freeTrajectoryState();
0122
0123   const GlobalVector& jetDirection(direction);
0124
0125   //
0126   // Check whether the track has been used in the vertex
0127   //
0128
0129   //FIXME
0130   float weight = 0.;  //vertex.trackWeight(aRecTrack);
0131
0132   TrajectoryStateOnSurface stateAtOrigin = track.impactPointState();
0133   TrajectoryStateOnSurface aTSOS = closestApproachToJet(*FTS, vertex, jetDirection, track.field());
0134   bool isValid = stateAtOrigin.isValid();
0135   //  bool IsValid= aTSOS.isValid();
0136
0137   if (isValid) {
0138     //get the Track line at origin
0139     Line::PositionType pos(stateAtOrigin.globalPosition());
0140     Line::DirectionType dir((stateAtOrigin.globalMomentum()).unit());
0141     Line track(pos, dir);
0142     // get the Jet  line
0143     // Vertex vertex(vertex);
0144     GlobalVector jetVector = jetDirection.unit();
0145     Line::PositionType pos2(GlobalPoint(vertex.x(), vertex.y(), vertex.z()));
0146     Line::DirectionType dir2(jetVector);
0147     Line jet(pos2, dir2);
0148     // now compute the distance between the two lines
0149     // If the track has been used to refit the Primary vertex then sign it positively, otherwise negative
0150
0151     theDistanceToJetAxis = (jet.distance(track)).mag();
0152     if (weight < 1)
0153       theDistanceToJetAxis = -theDistanceToJetAxis;
0154
0155     // ... and the flight distance along the Jet axis.
0156     GlobalPoint V = jet.position();
0157     GlobalVector Q = dir - jetVector.dot(dir) * jetVector;
0158     GlobalVector P = jetVector - jetVector.dot(dir) * dir;
0159     theDistanceAlongJetAxis = P.dot(V - pos) / Q.dot(dir);
0160
0161     //
0162     // get the covariance matrix of the vertex and compute the error on theDistanceToJetAxis
0163     //
0164
0165     ////AlgebraicSymMatrix vertexError = vertex.positionError().matrix();
0166
0167     // build the vector of closest approach between lines
0168
0169     GlobalVector H((jetVector.cross(dir).unit()));
0170
0171     CLHEP::HepVector Hh(3);
0172     Hh[0] = H.x();
0173     Hh[1] = H.y();
0174     Hh[2] = H.z();
0175
0176     //  theLDist_err = sqrt(vertexError.similarity(Hh));
0177
0178     //    cout << "distance to jet axis : "<< theDistanceToJetAxis <<" and error : "<< theLDist_err<<endl;
0179     // Now the impact parameter ...
0180
0181     /*    GlobalPoint T0 = track.position();
0182     GlobalVector D = (T0-V)- (T0-V).dot(dir) * dir;
0183     double IP = D.mag();
0184     GlobalVector Dold = distance(aTSOS, aJet.vertex(), jetDirection);
0185     double IPold = Dold.mag();
0186 */
0187   }
0188   Measurement1D DTJA(theDistanceToJetAxis, theLDist_err);
0189
0190   return pair<double, Measurement1D>(theDistanceAlongJetAxis, DTJA);
0191 }