SignedDecayLength3D.cc

Go to the documentation of this file.

#include <string>
 
#include "RecoBTag/BTagTools/interface/SignedDecayLength3D.h"
 
#include "DataFormats/GeometrySurface/interface/Line.h"
 
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/GeomPropagators/interface/AnalyticalTrajectoryExtrapolatorToLine.h"
#include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
 
#include "CLHEP/Vector/ThreeVector.h"
#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/Matrix/Vector.h"
 
using namespace std;
using namespace reco;
 
pair<bool, Measurement1D> SignedDecayLength3D::apply(const TransientTrack& transientTrack,
                                                     const GlobalVector& direction,
                                                     const Vertex& vertex) {
  double theError = 0.;
  bool theIsValid;
 
  //TrajectoryStateOnSurface TSOS = (aRecTrack).impactPointStateOnSurface();
  TrajectoryStateOnSurface TSOS = transientTrack.impactPointState();
  const FreeTrajectoryState* FTS = TSOS.freeTrajectoryState();
 
  TrajectoryStateOnSurface theTSOS = closestApproachToJet(*FTS, vertex, direction, transientTrack.field());
  theIsValid = theTSOS.isValid();
 
  if (theIsValid) {
    GlobalVector J = direction.unit();
    GlobalPoint vertexPosition(vertex.x(), vertex.y(), vertex.z());
 
    double theValue = J.dot(theTSOS.globalPosition() - vertexPosition);
 
    //error calculation
 
    AlgebraicVector3 j;
    j[0] = J.x();
    j[1] = J.y();
    j[2] = J.z();
    AlgebraicVector6 jj;
    jj[0] = J.x();
    jj[1] = J.y();
    jj[2] = J.z();
    jj[3] = 0.;
    jj[4] = 0.;
    jj[5] = 0.;  
    double E1 = ROOT::Math::Similarity(jj, theTSOS.cartesianError().matrix());
    //  double E2 = (aJet.vertex().positionError().matrix()).similarity(j);
    double E2 = ROOT::Math::Similarity(j, vertex.covariance());
 
    theError = sqrt(E1 + E2);
 
    //cout<< "Error ="<< theError<<endl;
    Measurement1D A(theValue, theError);
    return pair<bool, Measurement1D>(theIsValid, A);
  } else {
    return pair<bool, Measurement1D>(theIsValid, Measurement1D(0., 0.));
  }  // endif (isValid)
}  // end constructor declaration
 
TrajectoryStateOnSurface SignedDecayLength3D::closestApproachToJet(const FreeTrajectoryState& aFTS,
                                                                   const Vertex& vertex,
                                                                   const GlobalVector& aJetDirection,
                                                                   const MagneticField* field) {
  GlobalVector J = aJetDirection.unit();
 
  Line::PositionType pos(GlobalPoint(vertex.x(), vertex.y(), vertex.z()));
  Line::DirectionType dir(J);
  Line Jet(pos, dir);
 
  AnalyticalTrajectoryExtrapolatorToLine TETL(field);
 
  return TETL.extrapolate(aFTS, Jet);
}