PixelTrackBuilder.cc

Go to the documentation of this file.

#include "RecoPixelVertexing/PixelTrackFitting/interface/PixelTrackBuilder.h"


#include "DataFormats/GeometrySurface/interface/LocalError.h"
#include "DataFormats/GeometrySurface/interface/BoundPlane.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"

#include "TrackingTools/TrajectoryParametrization/interface/LocalTrajectoryError.h"
#include "TrackingTools/TrajectoryParametrization/interface/LocalTrajectoryParameters.h"
#include "TrackingTools/TrajectoryState/interface/BasicTrajectoryStateOnSurface.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"

#include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
#include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/TrackReco/interface/TrackBase.h"

#include <sstream>
using namespace std;
using namespace reco;


template <class T> T inline sqr( T t) {return t*t;}

namespace {
  __attribute__((unused))
  std::string print(
            const Measurement1D & pt,
            const Measurement1D & phi,
            const Measurement1D & cotTheta,
            const Measurement1D & tip,
            const Measurement1D & zip,
            float chi2,
            int   charge)
  {
    ostringstream str;
    str <<"\t pt: "  << pt.value() <<"+/-"<< pt.error()
        <<"\t phi: " << phi.value() <<"+/-"<< phi.error()
        <<"\t cot: " << cotTheta.value() <<"+/-"<< cotTheta.error()
        <<"\t tip: " << tip.value() <<"+/-"<< tip.error()
        <<"\t zip: " << zip.value() <<"+/-"<< zip.error()
        <<"\t chi2: " << chi2
        <<"\t charge: " << charge;
    return str.str();
  }

  __attribute__((unused))
  std::string print(const reco::Track & track, const GlobalPoint & origin)
  {

    math::XYZPoint bs(origin.x(), origin.y(), origin.z());

    Measurement1D phi( track.phi(), track.phiError());

    float theta = track.theta();
    float cosTheta = cos(theta);
    float sinTheta = sin(theta);
    float errLambda2 = sqr( track.lambdaError() );
    Measurement1D cotTheta(cosTheta/sinTheta, sqrt(errLambda2)/sqr(sinTheta));

    float pt_v = track.pt();
    float errInvP2 = sqr(track.qoverpError());
    float covIPtTheta = track.covariance(TrackBase::i_qoverp, TrackBase::i_lambda);
    float errInvPt2 = (   errInvP2
              + sqr(cosTheta/pt_v)*errLambda2
              + 2*(cosTheta/pt_v)*covIPtTheta
              ) / sqr(sinTheta);
    Measurement1D pt(pt_v, sqr(pt_v)*sqrt(errInvPt2));

    Measurement1D tip(track.dxy(bs), track.d0Error());

    Measurement1D zip(track.dz(bs), track.dzError());

    return print(pt, phi, cotTheta, tip, zip, track.chi2(), track.charge());
  }

  __attribute__((unused))
  std::string print(const BasicTrajectoryStateOnSurface & state)
  {
    // TrajectoryStateOnSurface state(bstate);
    float pt_v = state.globalMomentum().perp();
    float phi_v = state.globalMomentum().phi();
    float theta_v = state.globalMomentum().theta();

    CurvilinearTrajectoryError curv = state.curvilinearError();
    float errPhi2 = curv.matrix()(3,3);
    float errLambda2 = curv.matrix()(2,2);
    float errInvP2 = curv.matrix()(1,1);
    float covIPtTheta = curv.matrix()(1,2);
    float cosTheta = cos(theta_v);
    float sinTheta = sin(theta_v);
    float errInvPt2 = (   errInvP2
              + sqr(cosTheta/pt_v)*errLambda2
              + 2*(cosTheta/pt_v)*covIPtTheta) / sqr(sinTheta);
    float errCotTheta = sqrt(errLambda2)/sqr(sinTheta) ;
    Measurement1D pt(pt_v, sqr(pt_v) * sqrt(errInvPt2));
    Measurement1D phi(phi_v, sqrt(errPhi2) );
    Measurement1D cotTheta(cosTheta/sinTheta, errCotTheta);

    float zip_v = state.globalPosition().z();
    float zip_e = sqrt( state.localError().matrix()(4,4));
    Measurement1D zip(zip_v, zip_e);

    float tip_v  = state.localPosition().x();
    int tip_sign = (state.localMomentum().y()*cotTheta.value() > 0) ? -1 : 1;
    float tip_e  = sqrt( state.localError().matrix()(3,3) );
    Measurement1D tip( tip_sign*tip_v, tip_e);

    return print(pt, phi, cotTheta, tip, zip, 0., state.charge());
  }


#ifdef DEBUG_STATE
  inline void checkState(const  BasicTrajectoryStateOnSurface & bstate, const MagneticField* mf, const GlobalPoint & origin)
  {
    TrajectoryStateOnSurface state(bstate.clone());

    LogTrace("") << " *** PixelTrackBuilder::checkState: ";
    LogTrace("") << "INPUT,  ROTATION" << '\n' << state.surface().rotation();
    LogTrace("") << "INPUT,  TSOS:"<< '\n' << state;

    TransverseImpactPointExtrapolator tipe(mf);
    TrajectoryStateOnSurface test= tipe.extrapolate(state, origin);
    LogTrace("") << "CHECK-1 ROTATION" << '\n' << test.surface().rotation();
    LogTrace("") << "CHECK-1 TSOS" << '\n' << test;

    TSCPBuilderNoMaterial tscpBuilder;
    TrajectoryStateClosestToPoint tscp =
      tscpBuilder(*(state.freeState()), origin);
    const FreeTrajectoryState& fs = tscp.theState();
    LogTrace("") << "CHECK-2 FTS: " << fs;
  }
#endif

}


reco::Track * PixelTrackBuilder::build(
      const Measurement1D & pt,
      const Measurement1D & phi,
      const Measurement1D & cotTheta,
      const Measurement1D & tip,
      const Measurement1D & zip,
      float chi2,
      int   charge,
      const std::vector<const TrackingRecHit* >& hits,
      const MagneticField * mf,
      const GlobalPoint   & origin) const
{

  LogDebug("PixelTrackBuilder::build");
  LogTrace("") << "Reconstructed triplet kinematics: " << print(pt,phi,cotTheta,tip,zip,chi2,charge);

  double sinTheta = 1/std::sqrt(1+sqr(cotTheta.value()));
  double cosTheta = cotTheta.value()*sinTheta;
  int tipSign = tip.value() > 0 ? 1 : -1;

  AlgebraicSymMatrix55 m;
  double invPtErr = 1./sqr(pt.value()) * pt.error();
  m(0,0) = sqr(sinTheta) * (
              sqr(invPtErr)
            + sqr(cotTheta.error()/pt.value()*cosTheta * sinTheta)
            );
  m(0,2) = sqr( cotTheta.error()) * cosTheta * sqr(sinTheta) / pt.value();
  m(1,1) = sqr( phi.error() );
  m(2,2) = sqr( cotTheta.error());
  m(3,3) = sqr( tip.error() );
  m(4,4) = sqr( zip.error() );
  LocalTrajectoryError error(m);

  LocalTrajectoryParameters lpar(
    LocalPoint(tipSign*tip.value(), -tipSign*zip.value(), 0),
    LocalVector(0., -tipSign*pt.value()*cotTheta.value(), pt.value()),
    charge);


  float sp = std::sin(phi.value());
  float cp = std::cos(phi.value());
  Surface::RotationType rot(
                sp*tipSign, -cp*tipSign,           0,
                0         ,           0,    -tipSign,
                cp        ,  sp        ,           0);

  // BTSOS hold Surface in a shared pointer and  will be autodeleted when BTSOS goes out of scope...
  // to avoid memory churn we allocate it locally and just avoid it be deleted by refcount...
  Plane impPointPlane(origin, rot);
  // (twice just to be sure!)
  impPointPlane.addReference(); impPointPlane.addReference();
  // use Base (to avoid a useless new)
  BasicTrajectoryStateOnSurface impactPointState( lpar , error, impPointPlane, mf);

#ifdef DEBUG_STATE
  checkState(impactPointState,mf);
#endif
  LogTrace("") << "constructed TSOS:\n" << print(impactPointState);

  int ndof = 2*hits.size()-5;
  GlobalPoint vv = impactPointState.globalPosition();
  math::XYZPoint  pos( vv.x(), vv.y(), vv.z() );
  GlobalVector pp = impactPointState.globalMomentum();
  math::XYZVector mom( pp.x(), pp.y(), pp.z() );

  reco::Track * track = new reco::Track( chi2, ndof, pos, mom,
        impactPointState.charge(), impactPointState.curvilinearError());


  return track;
}