ParticleBuilder.cc

Go to the documentation of this file.

/* From SimpleFits Package
 * Designed an written by
 * author: Ian M. Nugent
 * Humboldt Foundations
 */
#include "RecoTauTag/ImpactParameter/interface/ErrorMatrixPropagator.h"
#include "RecoTauTag/ImpactParameter/interface/TrackHelixVertexFitter.h"
#include "RecoTauTag/ImpactParameter/interface/ParticleBuilder.h"
#include "RecoTauTag/ImpactParameter/interface/PDGInfo.h"
#include "Validation/EventGenerator/interface/PdtPdgMini.h"
#include "TrackingTools/TrajectoryParametrization/interface/PerigeeTrajectoryParameters.h"
#include "TrackingTools/TrajectoryParametrization/interface/PerigeeTrajectoryError.h"
#include "RecoTauTag/ImpactParameter/interface/TrackHelixVertexFitter.h"
#include <TVector3.h>
#include "FWCore/MessageLogger/interface/MessageLogger.h"

using namespace tauImpactParameter;

LorentzVectorParticle ParticleBuilder::createLorentzVectorParticle(const reco::TransientTrack& transTrk,
                                                                   const reco::Vertex& V,
                                                                   bool fromPerigee,
                                                                   bool useTrackHelixPropagation) {
  GlobalPoint p(V.position().x(), V.position().y(), V.position().z());
  TrackParticle tp = createTrackParticle(transTrk, p, fromPerigee, useTrackHelixPropagation);

  int N =
      TrackHelixVertexFitter::NFreeTrackPar + TrackHelixVertexFitter::NExtraPar + TrackHelixVertexFitter::MassOffSet;
  TVectorT<double> par(N);
  par(TrackHelixVertexFitter::x0) = V.position().x();
  par(TrackHelixVertexFitter::y0) = V.position().y();
  par(TrackHelixVertexFitter::z0) = V.position().z();
  par(TrackHelixVertexFitter::kappa0) = tp.parameter(TrackParticle::kappa);
  par(TrackHelixVertexFitter::lambda0) = tp.parameter(TrackParticle::lambda);
  par(TrackHelixVertexFitter::phi0) = tp.parameter(TrackParticle::phi);
  par(TrackHelixVertexFitter::NFreeTrackPar + TrackHelixVertexFitter::MassOffSet) = tp.mass();
  par(TrackHelixVertexFitter::NFreeTrackPar + TrackHelixVertexFitter::BField0) = tp.bField();

  TMatrixTSym<double> parcov(N);
  for (int i = 0; i < TrackHelixVertexFitter::NFreeVertexPar; i++) {
    for (int j = 0; j < TrackHelixVertexFitter::NFreeVertexPar; j++) {
      parcov(i, j) = V.covariance(i, j);
    }
  }
  parcov(TrackHelixVertexFitter::kappa0, TrackHelixVertexFitter::kappa0) =
      tp.covariance(TrackParticle::kappa, TrackParticle::kappa);
  parcov(TrackHelixVertexFitter::lambda0, TrackHelixVertexFitter::lambda0) =
      tp.covariance(TrackParticle::lambda, TrackParticle::lambda);
  parcov(TrackHelixVertexFitter::phi0, TrackHelixVertexFitter::phi0) =
      tp.covariance(TrackParticle::phi, TrackParticle::phi);

  TVectorT<double> LVPar = TrackHelixVertexFitter::computeLorentzVectorPar(par);
  TMatrixTSym<double> LVCov =
      ErrorMatrixPropagator::propagateError(&TrackHelixVertexFitter::computeLorentzVectorPar, par, parcov);
  return LorentzVectorParticle(LVPar, LVCov, tp.pdgId(), tp.charge(), tp.bField());
}

TrackParticle ParticleBuilder::createTrackParticle(const reco::TransientTrack& transTrk,
                                                   const GlobalPoint& p,
                                                   bool fromPerigee,
                                                   bool useTrackHelixPropagation) {
  // Configured for CMSSW Tracks only
  TVectorT<double> par(TrackParticle::NHelixPar + 1);
  TMatrixTSym<double> cov(TrackParticle::NHelixPar + 1);
  TVectorT<double> SFpar(TrackParticle::NHelixPar);
  TMatrixTSym<double> SFcov(TrackParticle::NHelixPar);
  if (!fromPerigee) {
    for (int i = 0; i < TrackParticle::NHelixPar; i++) {
      par(i) = transTrk.track().parameter(i);
      for (int j = 0; j < TrackParticle::NHelixPar; j++) {
        cov(i, j) = transTrk.track().covariance(i, j);
      }
    }
    par(TrackParticle::NHelixPar) = transTrk.field()->inInverseGeV(p).z();
    SFpar = convertCMSSWTrackParToSFTrackPar(par);
    SFcov = ErrorMatrixPropagator::propagateError(&ParticleBuilder::convertCMSSWTrackParToSFTrackPar, par, cov);
  } else {
    GlobalPoint TrackIPPos = transTrk.impactPointTSCP().position();
    //GlobalPoint TrackIPOrigin=transTrk.impactPointTSCP().referencePoint();
    GlobalPoint origin(0.0, 0.0, 0.0);
    for (int i = 0; i < TrackParticle::NHelixPar; i++) {
      par(i) = transTrk.trajectoryStateClosestToPoint(origin).perigeeParameters().vector()(i);
      for (int j = 0; j < TrackParticle::NHelixPar; j++) {
        cov(i, j) = transTrk.trajectoryStateClosestToPoint(origin).perigeeError().covarianceMatrix()(i, j);
      }
    }
    par(TrackParticle::NHelixPar) = transTrk.field()->inInverseGeV(p).z();
    SFpar = convertCMSSWTrackPerigeeToSFTrackPar(par);
    SFcov = ErrorMatrixPropagator::propagateError(&ParticleBuilder::convertCMSSWTrackPerigeeToSFTrackPar, par, cov);
    if (useTrackHelixPropagation) {
      // correct dxy dz neglecting material and radiative corrections

      LogDebug("RecoTauImpactParameterParticleBuilder")
          << "Offical CMS dxy - " << par(TrackParticle::dxy) << " dz " << par(TrackParticle::dz) << " kappa "
          << par(reco::TrackBase::i_qoverp);
      LogDebug("RecoTauImpactParameterParticleBuilder")
          << "Offical CMS dxy - SimpleFits Format" << SFpar(TrackParticle::dxy) << " dz " << SFpar(TrackParticle::dz)
          << " kappa " << SFpar(reco::TrackBase::i_qoverp);

      double x, y, z, dxy, dz, s, kappa, lambda, phi;
      TVectorT<double> freehelix(TrackHelixVertexFitter::NFreeTrackPar);
      freehelix(TrackHelixVertexFitter::x0) = TrackIPPos.x();
      freehelix(TrackHelixVertexFitter::y0) = TrackIPPos.y();
      freehelix(TrackHelixVertexFitter::z0) = TrackIPPos.z();
      freehelix(TrackHelixVertexFitter::kappa0) = SFpar(TrackParticle::kappa);
      freehelix(TrackHelixVertexFitter::lambda0) = SFpar(TrackParticle::lambda);
      freehelix(TrackHelixVertexFitter::phi0) = SFpar(TrackParticle::phi);
      TrackHelixVertexFitter::computedxydz(freehelix, 0, kappa, lambda, phi, x, y, z, s, dxy, dz);
      SFpar(TrackParticle::dxy) = dxy;
      SFpar(TrackParticle::dz) = dz;
      LogDebug("RecoTauImpactParameterParticleBuilder") << "Found values dxy " << dxy << " dz " << dz;
      //exit(0);
    }
  }

  PDGInfo pdgInfo;
  double c = transTrk.charge();
  return TrackParticle(
      SFpar, SFcov, abs(PdtPdgMini::pi_plus) * c, pdgInfo.pi_mass(), c, transTrk.field()->inInverseGeV(p).z());
}

reco::Vertex ParticleBuilder::getVertex(const LorentzVectorParticle& p) {
  TVector3 v = p.vertex();
  TMatrixTSym<double> vcov = p.vertexCov();
  reco::Vertex::Point vp(v.X(), v.Y(), v.Z());
  reco::Vertex::Error ve;
  for (int i = 0; i < vcov.GetNrows(); i++) {
    for (int j = 0; j < vcov.GetNrows(); j++) {
      ve(i, j) = vcov(i, j);
    }
  }
  return reco::Vertex(vp, ve);
}

TVectorT<double> ParticleBuilder::convertCMSSWTrackParToSFTrackPar(const TVectorT<double>& inpar) {
  TVectorT<double> par(TrackParticle::NHelixPar);
  par(TrackParticle::kappa) = -1.0 * inpar(TrackParticle::NHelixPar) * inpar(reco::TrackBase::i_qoverp) /
                              fabs(cos(inpar(reco::TrackBase::i_lambda)));
  par(TrackParticle::lambda) = inpar(reco::TrackBase::i_lambda);
  par(TrackParticle::phi) = inpar(reco::TrackBase::i_phi);
  par(TrackParticle::dz) = inpar(reco::TrackBase::i_dsz) / fabs(cos(inpar(reco::TrackBase::i_lambda)));
  par(TrackParticle::dxy) = inpar(reco::TrackBase::i_dxy);
  return par;
}

TVectorT<double> ParticleBuilder::convertCMSSWTrackPerigeeToSFTrackPar(const TVectorT<double>& inpar) {
  TVectorT<double> par(TrackParticle::NHelixPar);
  par(TrackParticle::kappa) = inpar(aCurv);
  par(TrackParticle::lambda) = TMath::Pi() / 2 - inpar(aTheta);
  par(TrackParticle::phi) = inpar(aPhi);
  par(TrackParticle::dxy) = -inpar(aTip);
  par(TrackParticle::dz) = inpar(aLip);
  return par;
}