File indexing completed on 2023-03-17 11:16:59
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
0068
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 }
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());
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
0127
0128
0129
0130 float weight = 0.;
0131
0132 TrajectoryStateOnSurface stateAtOrigin = track.impactPointState();
0133 TrajectoryStateOnSurface aTSOS = closestApproachToJet(*FTS, vertex, jetDirection, track.field());
0134 bool isValid = stateAtOrigin.isValid();
0135
0136
0137 if (isValid) {
0138
0139 Line::PositionType pos(stateAtOrigin.globalPosition());
0140 Line::DirectionType dir((stateAtOrigin.globalMomentum()).unit());
0141 Line track(pos, dir);
0142
0143
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
0149
0150
0151 theDistanceToJetAxis = (jet.distance(track)).mag();
0152 if (weight < 1)
0153 theDistanceToJetAxis = -theDistanceToJetAxis;
0154
0155
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
0163
0164
0165
0166
0167
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
0177
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187 }
0188 Measurement1D DTJA(theDistanceToJetAxis, theLDist_err);
0189
0190 return pair<double, Measurement1D>(theDistanceAlongJetAxis, DTJA);
0191 }