BTagTools/src/SignedTransverseImpactParameter.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/SignedTransverseImpactParameter.h"
0007 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0008 #include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
0009
0010 #include "CLHEP/Vector/ThreeVector.h"
0011 #include "CLHEP/Vector/LorentzVector.h"
0012 #include "CLHEP/Matrix/Vector.h"
0013 #include <string>
0014
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016
0017 using namespace std;
0018 using namespace reco;
0019
0020 pair<bool, Measurement1D> SignedTransverseImpactParameter::apply(const TransientTrack& track,
0021                                                                  const GlobalVector& direction,
0022                                                                  const Vertex& vertex) const {
0023   TrajectoryStateOnSurface stateOnSurface = track.impactPointState();
0024
0025   const FreeTrajectoryState* FTS = stateOnSurface.freeState();  //aRecTrack.stateAtFirstPoint().freeTrajectoryState());
0026
0027   GlobalPoint vertexPosition(vertex.x(), vertex.y(), vertex.z());
0028
0029   double theValue = 0.;
0030   double theError = 0.;
0031   TransverseImpactPointExtrapolator TIPE(track.field());
0032   TrajectoryStateOnSurface TSOS = TIPE.extrapolate(*FTS, vertexPosition);
0033
0034   if (!TSOS.isValid()) {
0035     theValue = 0.;
0036     theError = 0.;
0037   } else {
0038     GlobalPoint D0(TSOS.globalPosition());
0039
0040     GlobalVector DD0(D0.x() - vertex.x(), D0.y() - vertex.y(), 0.);
0041     const GlobalVector& JetDir(direction);
0042     double ps = DD0.dot(JetDir);
0043     theValue = DD0.mag() * (ps / abs(ps));
0044
0045     //error calculation
0046
0047     AlgebraicVector6 deriv;
0048     AlgebraicVector3 deriv_v;
0049     GlobalVector dd0 = DD0.unit();  //check
0050
0051     deriv_v[0] = -dd0.x();
0052     deriv_v[1] = -dd0.y();
0053     deriv_v[2] = -dd0.z();
0054
0055     deriv[0] = dd0.x();
0056     deriv[1] = dd0.y();
0057     deriv[2] = dd0.z();
0058     deriv[3] = 0.;
0059     deriv[4] = 0.;
0060     deriv[5] = 0.;
0061     //cout << TSOS.cartesianError().matrix() << endl;
0062     //cout << deriv << endl;
0063     double E1 = ROOT::Math::Similarity(deriv, TSOS.cartesianError().matrix());
0064     double E2 = ROOT::Math::Similarity(deriv_v, vertex.covariance());
0065     // (aJet.vertex().positionError().matrix()).similarity(deriv_v);
0066     theError = sqrt(E1 + E2);
0067     LogDebug("BTagTools") << "Tk error : " << E1 << " , Vtx error : " << E2 << "  tot : " << theError;
0068
0069   }  //end if
0070
0071   bool x = true;
0072
0073   Measurement1D A(theValue, theError);
0074   return pair<bool, Measurement1D>(x, A);
0075 }  // end constructor declaration
0076
0077 pair<bool, Measurement1D> SignedTransverseImpactParameter::zImpactParameter(const TransientTrack& track,
0078                                                                             const GlobalVector& direction,
0079                                                                             const Vertex& vertex) const {
0080   TransverseImpactPointExtrapolator TIPE(track.field());
0081   TrajectoryStateOnSurface TSOS = track.impactPointState();
0082
0083   if (!TSOS.isValid()) {
0084     LogDebug("BTagTools") << "====>>>> SignedTransverseImpactParameter::zImpactParameter : TSOS not valid";
0085     return pair<bool, Measurement1D>(false, Measurement1D(0.0, 0.0));
0086   }
0087
0088   GlobalPoint PV(vertex.x(), vertex.y(), vertex.z());
0089   TrajectoryStateOnSurface statePCA = TIPE.extrapolate(TSOS, PV);
0090
0091   if (!statePCA.isValid()) {
0092     LogDebug("BTagTools") << "====>>>> SignedTransverseImpactParameter::zImpactParameter : statePCA not valid";
0093     return pair<bool, Measurement1D>(false, Measurement1D(0.0, 0.0));
0094   }
0095
0096   GlobalPoint PCA = statePCA.globalPosition();
0097
0098   // sign as in rphi
0099   GlobalVector PVPCA(PCA.x() - PV.x(), PCA.y() - PV.y(), 0.);
0100   const GlobalVector& JetDir(direction);
0101   double sign = PVPCA.dot(JetDir);
0102   sign /= fabs(sign);
0103
0104   // z IP
0105   double deltaZ = fabs(PCA.z() - PV.z()) * sign;
0106
0107   // error
0108   double errPvZ2 = vertex.covariance()(2, 2);
0109   //CW  cout << "CW number or rows and columns : " << statePCA.cartesianError().matrix().num_row() << " , "
0110   //CW                                             << statePCA.cartesianError().matrix().num_col() << endl ;
0111   double errTrackZ2 = statePCA.cartesianError().matrix()(2, 2);
0112   double errZ = sqrt(errPvZ2 + errTrackZ2);
0113
0114   // CW alt. -> gives the same!!
0115   //  double errTZAlt = statePCA.localError().matrix()[4][4] ;
0116   //CW
0117
0118   return pair<bool, Measurement1D>(true, Measurement1D(deltaZ, errZ));
0119 }