IPTools/src/IPTools.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 "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h"
0006 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0007 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0008 #include "TrackingTools/IPTools/interface/IPTools.h"
0009 #include "CLHEP/Vector/ThreeVector.h"
0010 #include "CLHEP/Vector/LorentzVector.h"
0011 #include "CLHEP/Matrix/Vector.h"
0012 #include <string>
0013
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0016 #include "RecoVertex/VertexTools/interface/VertexDistanceXY.h"
0017 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
0018
0019 namespace IPTools {
0020   using namespace std;
0021   using namespace reco;
0022
0023   std::pair<bool, Measurement1D> absoluteImpactParameter(const TrajectoryStateOnSurface& tsos,
0024                                                          const reco::Vertex& vertex,
0025                                                          VertexDistance& distanceComputer) {
0026     if (!tsos.isValid()) {
0027       return pair<bool, Measurement1D>(false, Measurement1D(0., 0.));
0028     }
0029     GlobalPoint refPoint = tsos.globalPosition();
0030     GlobalError refPointErr = tsos.cartesianError().position();
0031     GlobalPoint vertexPosition = RecoVertex::convertPos(vertex.position());
0032     GlobalError vertexPositionErr = RecoVertex::convertError(vertex.error());
0033     return pair<bool, Measurement1D>(
0034         true,
0035         distanceComputer.distance(VertexState(vertexPosition, vertexPositionErr), VertexState(refPoint, refPointErr)));
0036   }
0037
0038   std::pair<bool, Measurement1D> absoluteImpactParameter3D(const reco::TransientTrack& transientTrack,
0039                                                            const reco::Vertex& vertex) {
0040     AnalyticalImpactPointExtrapolator extrapolator(transientTrack.field());
0041     VertexDistance3D dist;
0042     return absoluteImpactParameter(
0043         extrapolator.extrapolate(transientTrack.impactPointState(), RecoVertex::convertPos(vertex.position())),
0044         vertex,
0045         dist);
0046   }
0047   std::pair<bool, Measurement1D> absoluteTransverseImpactParameter(const reco::TransientTrack& transientTrack,
0048                                                                    const reco::Vertex& vertex) {
0049     TransverseImpactPointExtrapolator extrapolator(transientTrack.field());
0050     VertexDistanceXY dist;
0051     return absoluteImpactParameter(
0052         extrapolator.extrapolate(transientTrack.impactPointState(), RecoVertex::convertPos(vertex.position())),
0053         vertex,
0054         dist);
0055   }
0056
0057   pair<bool, Measurement1D> signedTransverseImpactParameter(const TransientTrack& track,
0058                                                             const GlobalVector& direction,
0059                                                             const Vertex& vertex) {
0060     //Extrapolate to closest point on transverse plane
0061     TransverseImpactPointExtrapolator extrapolator(track.field());
0062     TrajectoryStateOnSurface closestOnTransversePlaneState =
0063         extrapolator.extrapolate(track.impactPointState(), RecoVertex::convertPos(vertex.position()));
0064
0065     //Compute absolute value
0066     VertexDistanceXY dist;
0067     pair<bool, Measurement1D> result = absoluteImpactParameter(closestOnTransversePlaneState, vertex, dist);
0068     if (!result.first)
0069       return result;
0070
0071     //Compute Sign
0072     GlobalPoint impactPoint = closestOnTransversePlaneState.globalPosition();
0073     GlobalVector IPVec(impactPoint.x() - vertex.x(), impactPoint.y() - vertex.y(), 0.);
0074     double prod = IPVec.dot(direction);
0075     double sign = (prod >= 0) ? 1. : -1.;
0076
0077     //Apply sign to the result
0078     return pair<bool, Measurement1D>(result.first, Measurement1D(sign * result.second.value(), result.second.error()));
0079   }
0080
0081   pair<bool, Measurement1D> signedImpactParameter3D(const TransientTrack& track,
0082                                                     const GlobalVector& direction,
0083                                                     const Vertex& vertex) {
0084     //Extrapolate to closest point on transverse plane
0085     AnalyticalImpactPointExtrapolator extrapolator(track.field());
0086     TrajectoryStateOnSurface closestIn3DSpaceState =
0087         extrapolator.extrapolate(track.impactPointState(), RecoVertex::convertPos(vertex.position()));
0088
0089     //Compute absolute value
0090     VertexDistance3D dist;
0091     pair<bool, Measurement1D> result = absoluteImpactParameter(closestIn3DSpaceState, vertex, dist);
0092     if (!result.first)
0093       return result;
0094
0095     //Compute Sign
0096     GlobalPoint impactPoint = closestIn3DSpaceState.globalPosition();
0097     GlobalVector IPVec(impactPoint.x() - vertex.x(), impactPoint.y() - vertex.y(), impactPoint.z() - vertex.z());
0098     double prod = IPVec.dot(direction);
0099     double sign = (prod >= 0) ? 1. : -1.;
0100
0101     //Apply sign to the result
0102     return pair<bool, Measurement1D>(result.first, Measurement1D(sign * result.second.value(), result.second.error()));
0103   }
0104
0105   pair<bool, Measurement1D> signedDecayLength3D(const TrajectoryStateOnSurface& closestToJetState,
0106                                                 const GlobalVector& direction,
0107                                                 const Vertex& vertex) {
0108     //Check if extrapolation has been successfull
0109     if (!closestToJetState.isValid()) {
0110       return pair<bool, Measurement1D>(false, Measurement1D(0., 0.));
0111     }
0112
0113     GlobalVector jetDirection = direction.unit();
0114     GlobalPoint vertexPosition(vertex.x(), vertex.y(), vertex.z());
0115
0116     double decayLen = jetDirection.dot(closestToJetState.globalPosition() - vertexPosition);
0117
0118     //error calculation
0119
0120     AlgebraicVector3 j;
0121     j[0] = jetDirection.x();
0122     j[1] = jetDirection.y();
0123     j[2] = jetDirection.z();
0124     AlgebraicVector6 jj;
0125     jj[0] = jetDirection.x();
0126     jj[1] = jetDirection.y();
0127     jj[2] = jetDirection.z();
0128     jj[3] = 0.;
0129     jj[4] = 0.;
0130     jj[5] = 0.;  ///chech it!!!!!!!!!!!!!!!!!!!!!!!
0131
0132     //TODO: FIXME: the extrapolation uncertainty is very relevant here should be taken into account!!
0133     double trackError2 = ROOT::Math::Similarity(jj, closestToJetState.cartesianError().matrix());
0134     double vertexError2 = ROOT::Math::Similarity(j, vertex.covariance());
0135
0136     double decayLenError = sqrt(trackError2 + vertexError2);
0137
0138     return pair<bool, Measurement1D>(true, Measurement1D(decayLen, decayLenError));
0139   }
0140
0141   pair<bool, Measurement1D> linearizedSignedImpactParameter3D(const TrajectoryStateOnSurface& closestToJetState,
0142                                                               const GlobalVector& direction,
0143                                                               const Vertex& vertex) {
0144     //Check if extrapolation has been successfull
0145     if (!closestToJetState.isValid()) {
0146       return pair<bool, Measurement1D>(false, Measurement1D(0., 0.));
0147     }
0148
0149     GlobalPoint vertexPosition(vertex.x(), vertex.y(), vertex.z());
0150     GlobalVector impactParameter = linearImpactParameter(closestToJetState, vertexPosition);
0151     GlobalVector jetDir = direction.unit();
0152     GlobalVector flightDistance(closestToJetState.globalPosition() - vertexPosition);
0153     double theDistanceAlongJetAxis = jetDir.dot(flightDistance);
0154     double signedIP = impactParameter.mag() *
0155                       ((theDistanceAlongJetAxis != 0) ? theDistanceAlongJetAxis / abs(theDistanceAlongJetAxis) : 1.);
0156
0157     GlobalVector ipDirection = impactParameter.unit();
0158     //GlobalPoint closestPoint = closestToJetState.globalPosition();
0159     GlobalVector momentumAtClosestPoint = closestToJetState.globalMomentum();
0160     GlobalVector momentumDir = momentumAtClosestPoint.unit();
0161
0162     AlgebraicVector3 deriv_v;
0163     deriv_v[0] = -ipDirection.x();
0164     deriv_v[1] = -ipDirection.y();
0165     deriv_v[2] = -ipDirection.z();
0166
0167     AlgebraicVector6 deriv;
0168     deriv[0] = ipDirection.x();
0169     deriv[1] = ipDirection.y();
0170     deriv[2] = ipDirection.z();
0171     deriv[3] = -(momentumDir.dot(flightDistance) * ipDirection.x()) / momentumAtClosestPoint.mag();
0172     deriv[4] = -(momentumDir.dot(flightDistance) * ipDirection.y()) / momentumAtClosestPoint.mag();
0173     deriv[5] = -(momentumDir.dot(flightDistance) * ipDirection.z()) / momentumAtClosestPoint.mag();
0174
0175     double trackError2 = ROOT::Math::Similarity(deriv, closestToJetState.cartesianError().matrix());
0176     double vertexError2 = ROOT::Math::Similarity(deriv_v, vertex.covariance());
0177     double ipError = sqrt(trackError2 + vertexError2);
0178
0179     return pair<bool, Measurement1D>(true, Measurement1D(signedIP, ipError));
0180   }
0181
0182   TrajectoryStateOnSurface closestApproachToJet(const TrajectoryStateOnSurface& state,
0183                                                 const Vertex& vertex,
0184                                                 const GlobalVector& direction,
0185                                                 const MagneticField* field) {
0186     Line::PositionType pos(GlobalPoint(vertex.x(), vertex.y(), vertex.z()));
0187     Line::DirectionType dir(direction.unit());
0188     Line jetLine(pos, dir);
0189
0190     AnalyticalTrajectoryExtrapolatorToLine extrapolator(field);
0191
0192     return extrapolator.extrapolate(state, jetLine);
0193   }
0194
0195   /**
0196    * Compute the impact parameter of a track, linearized from the given state, with respect to a given point
0197    */
0198   GlobalVector linearImpactParameter(const TrajectoryStateOnSurface& state, const GlobalPoint& point) {
0199     Line::PositionType pos(state.globalPosition());
0200     Line::DirectionType dir((state.globalMomentum()).unit());
0201     Line trackLine(pos, dir);
0202     const GlobalPoint& tmp = point;
0203     return trackLine.distance(tmp);
0204   }
0205
0206   pair<double, Measurement1D> jetTrackDistance(const TransientTrack& track,
0207                                                const GlobalVector& direction,
0208                                                const Vertex& vertex) {
0209     double theLDist_err(0.);
0210
0211     //FIXME
0212     float weight = 0.;  //vertex.trackWeight(track);
0213
0214     TrajectoryStateOnSurface stateAtOrigin = track.impactPointState();
0215     if (!stateAtOrigin.isValid()) {
0216       //TODO: throw instead?
0217       return pair<bool, Measurement1D>(false, Measurement1D(0., 0.));
0218     }
0219
0220     //get the Track line at origin
0221     Line::PositionType posTrack(stateAtOrigin.globalPosition());
0222     Line::DirectionType dirTrack((stateAtOrigin.globalMomentum()).unit());
0223     Line trackLine(posTrack, dirTrack);
0224     // get the Jet  line
0225     // Vertex vertex(vertex);
0226     GlobalVector jetVector = direction.unit();
0227     Line::PositionType posJet(GlobalPoint(vertex.x(), vertex.y(), vertex.z()));
0228     Line::DirectionType dirJet(jetVector);
0229     Line jetLine(posJet, dirJet);
0230
0231     // now compute the distance between the two lines
0232     // If the track has been used to refit the Primary vertex then sign it positively, otherwise negative
0233     double theDistanceToJetAxis = (jetLine.distance(trackLine)).mag();
0234     if (weight < 1)
0235       theDistanceToJetAxis = -theDistanceToJetAxis;
0236
0237     // ... and the flight distance along the Jet axis.
0238     GlobalPoint V = jetLine.position();
0239     GlobalVector Q = dirTrack - jetVector.dot(dirTrack) * jetVector;
0240     GlobalVector P = jetVector - jetVector.dot(dirTrack) * dirTrack;
0241     double theDistanceAlongJetAxis = P.dot(V - posTrack) / Q.dot(dirTrack);
0242
0243     //
0244     // get the covariance matrix of the vertex and compute the error on theDistanceToJetAxis
0245     //
0246
0247     ////AlgebraicSymMatrix vertexError = vertex.positionError().matrix();
0248
0249     // build the vector of closest approach between lines
0250
0251     //FIXME: error not computed.
0252     GlobalVector H((jetVector.cross(dirTrack).unit()));
0253     CLHEP::HepVector Hh(3);
0254     Hh[0] = H.x();
0255     Hh[1] = H.y();
0256     Hh[2] = H.z();
0257
0258     //  theLDist_err = sqrt(vertexError.similarity(Hh));
0259
0260     //    cout << "distance to jet axis : "<< theDistanceToJetAxis <<" and error : "<< theLDist_err<<endl;
0261     // Now the impact parameter ...
0262
0263     /*    GlobalPoint T0 = track.position();
0264       GlobalVector D = (T0-V)- (T0-V).dot(dir) * dir;
0265       double IP = D.mag();
0266       GlobalVector Dold = distance(aTSOS, aJet.vertex(), jetDirection);
0267       double IPold = Dold.mag();
0268     */
0269
0270     Measurement1D DTJA(theDistanceToJetAxis, theLDist_err);
0271
0272     return pair<double, Measurement1D>(theDistanceAlongJetAxis, DTJA);
0273   }
0274
0275 }  // namespace IPTools