/data/refman/pasoursint/CMSSW_4_2_9_HLT1_bphpatch4/src/TrackingTools/IPTools/src/IPTools.cc

Go to the documentation of this file.

#include "DataFormats/GeometrySurface/interface/Line.h"

#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/GeomPropagators/interface/AnalyticalTrajectoryExtrapolatorToLine.h"
#include "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h"
#include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
#include "TrackingTools/IPTools/interface/IPTools.h"
#include "CLHEP/Vector/ThreeVector.h"
#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/Matrix/Vector.h"
#include <string>

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
#include "RecoVertex/VertexTools/interface/VertexDistanceXY.h"
#include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"


namespace IPTools
{
  using namespace std;
  using namespace reco;


    std::pair<bool,Measurement1D> absoluteImpactParameter(const TrajectoryStateOnSurface & tsos  , const  reco::Vertex & vertex, VertexDistance & distanceComputer) {
        if(!tsos.isValid()) {
         return pair<bool,Measurement1D>(false,Measurement1D(0.,0.)) ;
        }
        GlobalPoint refPoint          = tsos.globalPosition();
        GlobalError refPointErr       = tsos.cartesianError().position();
        GlobalPoint vertexPosition    = RecoVertex::convertPos(vertex.position());
        GlobalError vertexPositionErr = RecoVertex::convertError(vertex.error());
        return pair<bool,Measurement1D>(true,distanceComputer.distance(VertexState(vertexPosition,vertexPositionErr), VertexState(refPoint, refPointErr)));
   }

    std::pair<bool,Measurement1D> absoluteImpactParameter3D(const reco::TransientTrack & transientTrack, const  reco::Vertex & vertex)
    {
      AnalyticalImpactPointExtrapolator extrapolator(transientTrack.field());
      VertexDistance3D dist;
      return absoluteImpactParameter(extrapolator.extrapolate(transientTrack.impactPointState(), RecoVertex::convertPos(vertex.position())), vertex, dist);
    }
    std::pair<bool,Measurement1D> absoluteTransverseImpactParameter(const reco::TransientTrack & transientTrack, const  reco::Vertex & vertex)
    {
      TransverseImpactPointExtrapolator extrapolator(transientTrack.field());
      VertexDistanceXY dist;
      return absoluteImpactParameter(extrapolator.extrapolate(transientTrack.impactPointState(), RecoVertex::convertPos(vertex.position())), vertex,  dist);
    }
    
  pair<bool,Measurement1D> signedTransverseImpactParameter(const TransientTrack & track,
                                                           const GlobalVector & direction, const  Vertex & vertex){
    //Extrapolate to closest point on transverse plane
    TransverseImpactPointExtrapolator extrapolator(track.field());
    TrajectoryStateOnSurface closestOnTransversePlaneState = extrapolator.extrapolate(track.impactPointState(),RecoVertex::convertPos(vertex.position()));
    
    //Compute absolute value
    VertexDistanceXY dist;
    pair<bool,Measurement1D> result = absoluteImpactParameter(closestOnTransversePlaneState, vertex,  dist);
    if(!result.first) return result;

    //Compute Sign
    GlobalPoint impactPoint    = closestOnTransversePlaneState.globalPosition();
    GlobalVector IPVec(impactPoint.x()-vertex.x(),impactPoint.y()-vertex.y(),0.);
    double prod = IPVec.dot(direction);
    double sign = (prod>=0) ? 1. : -1.;
    
    //Apply sign to the result
    return pair<bool,Measurement1D>(result.first,Measurement1D(sign*result.second.value(), result.second.error()));
  }

  pair<bool,Measurement1D> signedImpactParameter3D(const TransientTrack & track,
                                                           const GlobalVector & direction, const  Vertex & vertex){
    //Extrapolate to closest point on transverse plane
    AnalyticalImpactPointExtrapolator extrapolator(track.field());
    TrajectoryStateOnSurface closestIn3DSpaceState = extrapolator.extrapolate(track.impactPointState(),RecoVertex::convertPos(vertex.position()));

    //Compute absolute value
    VertexDistance3D dist;
    pair<bool,Measurement1D> result = absoluteImpactParameter(closestIn3DSpaceState, vertex,  dist);
    if(!result.first) return result;

    //Compute Sign
    GlobalPoint impactPoint = closestIn3DSpaceState.globalPosition();
    GlobalVector IPVec(impactPoint.x()-vertex.x(),impactPoint.y()-vertex.y(),impactPoint.z()-vertex.z());
    double prod = IPVec.dot(direction);
    double sign = (prod>=0) ? 1. : -1.;

    //Apply sign to the result
    return pair<bool,Measurement1D>(result.first,Measurement1D(sign*result.second.value(), result.second.error()));
  }



  pair<bool,Measurement1D> signedDecayLength3D(const   TrajectoryStateOnSurface & closestToJetState,
                                               const GlobalVector & direction, const  Vertex & vertex)  {

    //Check if extrapolation has been successfull
    if(!closestToJetState.isValid()) {
      return pair<bool,Measurement1D>(false,Measurement1D(0.,0.));
    }

    GlobalVector jetDirection = direction.unit();
    GlobalPoint vertexPosition(vertex.x(),vertex.y(),vertex.z());
  
    double decayLen = jetDirection.dot(closestToJetState.globalPosition()-vertexPosition);

    //error calculation

    AlgebraicVector3 j;
    j[0] = jetDirection.x();
    j[1] = jetDirection.y();
    j[2] = jetDirection.z();
    AlgebraicVector6 jj;
    jj[0] = jetDirection.x();
    jj[1] = jetDirection.y();
    jj[2] = jetDirection.z();
    jj[3] =0.;
    jj[4] =0.;
    jj[5] =0.;
     
    //TODO: FIXME: the extrapolation uncertainty is very relevant here should be taken into account!!
    double trackError2 = ROOT::Math::Similarity(jj,closestToJetState.cartesianError().matrix());
    double vertexError2 = ROOT::Math::Similarity(j,vertex.covariance());

    double decayLenError = sqrt(trackError2+vertexError2);
    
    return pair<bool,Measurement1D>(true,Measurement1D(decayLen,decayLenError));

  }


    
  pair<bool,Measurement1D> linearizedSignedImpactParameter3D(const   TrajectoryStateOnSurface & closestToJetState ,
                                                   const GlobalVector & direction, const  Vertex & vertex)
  {
    //Check if extrapolation has been successfull
    if(!closestToJetState.isValid()) {
      return pair<bool,Measurement1D>(false,Measurement1D(0.,0.));
    }

    GlobalPoint vertexPosition(vertex.x(),vertex.y(),vertex.z());
    GlobalVector impactParameter = linearImpactParameter(closestToJetState, vertexPosition);
    GlobalVector jetDir = direction.unit();
    GlobalVector flightDistance(closestToJetState.globalPosition()-vertexPosition);
    double theDistanceAlongJetAxis = jetDir.dot(flightDistance);
    double signedIP = impactParameter.mag()*((theDistanceAlongJetAxis!=0)?theDistanceAlongJetAxis/abs(theDistanceAlongJetAxis):1.);


    GlobalVector ipDirection = impactParameter.unit();
    GlobalPoint closestPoint = closestToJetState.globalPosition();
    GlobalVector momentumAtClosestPoint = closestToJetState.globalMomentum(); 
    GlobalVector momentumDir = momentumAtClosestPoint.unit();
    
    AlgebraicVector3 deriv_v;
    deriv_v[0] = - ipDirection.x();
    deriv_v[1] = - ipDirection.y();
    deriv_v[2] = - ipDirection.z();

 
    AlgebraicVector6 deriv;
    deriv[0] = ipDirection.x();
    deriv[1] = ipDirection.y();
    deriv[2] = ipDirection.z();
    deriv[3] = -  (momentumDir.dot(flightDistance)*ipDirection.x())/momentumAtClosestPoint.mag();
    deriv[4] = -  (momentumDir.dot(flightDistance)*ipDirection.y())/momentumAtClosestPoint.mag();
    deriv[5] = -  (momentumDir.dot(flightDistance)*ipDirection.z())/momentumAtClosestPoint.mag();

    double trackError2 = ROOT::Math::Similarity(deriv , closestToJetState.cartesianError().matrix());
    double vertexError2 = ROOT::Math::Similarity(deriv_v , vertex.covariance());
    double ipError = sqrt(trackError2+vertexError2);

    return pair<bool,Measurement1D>(true,Measurement1D(signedIP,ipError));
  }



  TrajectoryStateOnSurface closestApproachToJet(const TrajectoryStateOnSurface & state,const Vertex & vertex, const GlobalVector& direction,const MagneticField * field) {
  
    Line::PositionType pos(GlobalPoint(vertex.x(),vertex.y(),vertex.z()));
    Line::DirectionType dir(direction.unit());
    Line jetLine(pos,dir);
  
    AnalyticalTrajectoryExtrapolatorToLine extrapolator(field);

    return extrapolator.extrapolate(state, jetLine);
  }

  GlobalVector linearImpactParameter(const TrajectoryStateOnSurface & state, const GlobalPoint & point)  {

    Line::PositionType pos(state.globalPosition());
    Line::DirectionType dir((state.globalMomentum()).unit());
    Line trackLine(pos,dir);
    GlobalPoint  tmp=point; 
    return  trackLine.distance(tmp);
  }

  pair<double,Measurement1D> jetTrackDistance(const TransientTrack & track, const GlobalVector & direction, const Vertex & vertex) {
    double  theLDist_err(0.);
  
    //FIXME
    float weight=0.;//vertex.trackWeight(track);

    TrajectoryStateOnSurface stateAtOrigin = track.impactPointState(); 
    if(!stateAtOrigin.isValid())
      {
        //TODO: throw instead?
        return pair<bool,Measurement1D>(false,Measurement1D(0.,0.));
      }
   
    //get the Track line at origin
    Line::PositionType posTrack(stateAtOrigin.globalPosition());
    Line::DirectionType dirTrack((stateAtOrigin.globalMomentum()).unit());
    Line trackLine(posTrack,dirTrack);
    // get the Jet  line 
    // Vertex vertex(vertex);
    GlobalVector jetVector = direction.unit();    
    Line::PositionType posJet(GlobalPoint(vertex.x(),vertex.y(),vertex.z()));
    Line::DirectionType dirJet(jetVector);
    Line jetLine(posJet,dirJet);
  
    // now compute the distance between the two lines
    // If the track has been used to refit the Primary vertex then sign it positively, otherwise negative
    double theDistanceToJetAxis = (jetLine.distance(trackLine)).mag();
    if (weight<1) theDistanceToJetAxis= -theDistanceToJetAxis;

    // ... and the flight distance along the Jet axis.
    GlobalPoint  V = jetLine.position();    
    GlobalVector Q = dirTrack - jetVector.dot(dirTrack) * jetVector;
    GlobalVector P = jetVector - jetVector.dot(dirTrack) * dirTrack;
    double theDistanceAlongJetAxis = P.dot(V-posTrack)/Q.dot(dirTrack);

    //
    // get the covariance matrix of the vertex and compute the error on theDistanceToJetAxis
    //
    

    // build the vector of closest approach between lines


    //FIXME: error not computed.
    GlobalVector H((jetVector.cross(dirTrack).unit()));
    CLHEP::HepVector Hh(3);
    Hh[0] = H.x();
    Hh[1] = H.y();
    Hh[2] = H.z();
    
    //  theLDist_err = sqrt(vertexError.similarity(Hh));

    //    cout << "distance to jet axis : "<< theDistanceToJetAxis <<" and error : "<< theLDist_err<<endl;
    // Now the impact parameter ...

    /*    GlobalPoint T0 = track.position();
          GlobalVector D = (T0-V)- (T0-V).dot(dir) * dir;
          double IP = D.mag();    
          GlobalVector Dold = distance(aTSOS, aJet.vertex(), jetDirection);
          double IPold = Dold.mag();
    */



  
    Measurement1D DTJA(theDistanceToJetAxis,theLDist_err);
  
    return pair<double,Measurement1D> (theDistanceAlongJetAxis,DTJA);
  }

}