BTagTools/src/SignedDecayLength3D.cc

0001 #include <string>
0002
0003 #include "RecoBTag/BTagTools/interface/SignedDecayLength3D.h"
0004
0005 #include "DataFormats/GeometrySurface/interface/Line.h"
0006
0007 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0008 #include "TrackingTools/GeomPropagators/interface/AnalyticalTrajectoryExtrapolatorToLine.h"
0009 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0010 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0011
0012 #include "CLHEP/Vector/ThreeVector.h"
0013 #include "CLHEP/Vector/LorentzVector.h"
0014 #include "CLHEP/Matrix/Vector.h"
0015
0016 using namespace std;
0017 using namespace reco;
0018
0019 pair<bool, Measurement1D> SignedDecayLength3D::apply(const TransientTrack& transientTrack,
0020                                                      const GlobalVector& direction,
0021                                                      const Vertex& vertex) {
0022   double theError = 0.;
0023   bool theIsValid;
0024
0025   //TrajectoryStateOnSurface TSOS = (aRecTrack).impactPointStateOnSurface();
0026   TrajectoryStateOnSurface TSOS = transientTrack.impactPointState();
0027   const FreeTrajectoryState* FTS = TSOS.freeTrajectoryState();
0028
0029   TrajectoryStateOnSurface theTSOS = closestApproachToJet(*FTS, vertex, direction, transientTrack.field());
0030   theIsValid = theTSOS.isValid();
0031
0032   if (theIsValid) {
0033     GlobalVector J = direction.unit();
0034     GlobalPoint vertexPosition(vertex.x(), vertex.y(), vertex.z());
0035
0036     double theValue = J.dot(theTSOS.globalPosition() - vertexPosition);
0037
0038     //error calculation
0039
0040     AlgebraicVector3 j;
0041     j[0] = J.x();
0042     j[1] = J.y();
0043     j[2] = J.z();
0044     AlgebraicVector6 jj;
0045     jj[0] = J.x();
0046     jj[1] = J.y();
0047     jj[2] = J.z();
0048     jj[3] = 0.;
0049     jj[4] = 0.;
0050     jj[5] = 0.;  ///chech it!!!!!!!!!!!!!!!!!!!!!!!
0051     double E1 = ROOT::Math::Similarity(jj, theTSOS.cartesianError().matrix());
0052     //  double E2 = (aJet.vertex().positionError().matrix()).similarity(j);
0053     double E2 = ROOT::Math::Similarity(j, vertex.covariance());
0054
0055     theError = sqrt(E1 + E2);
0056
0057     //cout<< "Error ="<< theError<<endl;
0058     Measurement1D A(theValue, theError);
0059     return pair<bool, Measurement1D>(theIsValid, A);
0060   } else {
0061     return pair<bool, Measurement1D>(theIsValid, Measurement1D(0., 0.));
0062   }  // endif (isValid)
0063 }  // end constructor declaration
0064
0065 TrajectoryStateOnSurface SignedDecayLength3D::closestApproachToJet(const FreeTrajectoryState& aFTS,
0066                                                                    const Vertex& vertex,
0067                                                                    const GlobalVector& aJetDirection,
0068                                                                    const MagneticField* field) {
0069   GlobalVector J = aJetDirection.unit();
0070
0071   Line::PositionType pos(GlobalPoint(vertex.x(), vertex.y(), vertex.z()));
0072   Line::DirectionType dir(J);
0073   Line Jet(pos, dir);
0074
0075   AnalyticalTrajectoryExtrapolatorToLine TETL(field);
0076
0077   return TETL.extrapolate(aFTS, Jet);
0078 }