TrackHelixVertexFitter.cc

Go to the documentation of this file.

/* From SimpleFits Package
 * Designed an written by
 * author: Ian M. Nugent
 * Humboldt Foundations
 */
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "RecoTauTag/ImpactParameter/interface/TrackHelixVertexFitter.h"
#include "TDecompBK.h"
#include <iostream>

using namespace tauImpactParameter;

TrackHelixVertexFitter::TrackHelixVertexFitter(const std::vector<TrackParticle>& particles, const TVector3& vguess)
    : isFit_(false),
      isConfigured_(false),
      nParticles_(particles.size()),
      nPar_((NFreeTrackPar - NFreeVertexPar) * particles.size() + NFreeVertexPar),
      nVal_(TrackParticle::NHelixPar * particles.size()) {
  particles_ = particles;
  par_.ResizeTo(nPar_);
  parcov_.ResizeTo(nPar_, nPar_);
  val_.ResizeTo(nVal_);
  cov_.ResizeTo(nVal_, nVal_);
  for (unsigned int p = 0; p < particles.size(); p++) {
    for (unsigned int j = 0; j < TrackParticle::NHelixPar; j++) {
      val_(measuredValueIndex(j, p)) = particles[p].parameter(j);
      for (unsigned int k = 0; k < TrackParticle::NHelixPar; k++) {
        cov_(measuredValueIndex(j, p), measuredValueIndex(k, p)) = particles[p].covariance(j, k);
      }
    }
  }
  TDecompBK Inverter(cov_);
  double det = cov_.Determinant();
  if (!Inverter.Decompose()) {
    edm::LogWarning("TrackHelixVertexFitter::TrackHelixVertexFitter")
        << "Fit failed: unable to invert SYM gain matrix " << det << " \n"
        << std::endl;
    return;
  }

  cov_inv_.ResizeTo(nVal_, nVal_);
  cov_inv_ = Inverter.Invert();
  ndf_ = nVal_ - nPar_;
  // Set Initial conditions within reason
  par_(x0) = vguess.X();
  parcov_(x0, x0) = 1.0;
  par_(y0) = vguess.Y();
  parcov_(y0, y0) = 1.0;
  par_(z0) = vguess.Z();
  parcov_(z0, z0) = 1.0;
  for (unsigned int p = 0; p < particles_.size(); p++) {
    par_(freeParIndex(kappa0, p)) = val_(measuredValueIndex(TrackParticle::kappa, p));
    par_(freeParIndex(lambda0, p)) = val_(measuredValueIndex(TrackParticle::lambda, p));
    par_(freeParIndex(phi0, p)) = val_(measuredValueIndex(TrackParticle::phi, p));

    parcov_(freeParIndex(kappa0, p), freeParIndex(kappa0, p)) =
        cov_(measuredValueIndex(TrackParticle::kappa, p), measuredValueIndex(TrackParticle::kappa, p));
    parcov_(freeParIndex(lambda0, p), freeParIndex(lambda0, p)) =
        cov_(measuredValueIndex(TrackParticle::lambda, p), measuredValueIndex(TrackParticle::lambda, p));
    parcov_(freeParIndex(phi0, p), freeParIndex(phi0, p)) =
        cov_(measuredValueIndex(TrackParticle::phi, p), measuredValueIndex(TrackParticle::phi, p));
  }
  isConfigured_ = true;
}

TrackHelixVertexFitter::~TrackHelixVertexFitter() {}

double TrackHelixVertexFitter::updateChisquare(const TVectorT<double>& inpar) {
  TVectorT<double> vprime = computePar(inpar);
  TVectorT<double> dalpha = vprime - val_;
  double c2 = dalpha * (cov_inv_ * dalpha);
  return c2;
}

std::vector<TrackParticle> TrackHelixVertexFitter::getRefitTracks() {
  std::vector<TrackParticle> refitParticles;
  for (unsigned int p = 0; p < particles_.size(); p++) {
    TVectorT<double> FreePar(NFreeTrackPar);
    TMatrixTSym<double> FreeParCov(NFreeTrackPar);
    for (int i = 0; i < FreeParCov.GetNrows(); i++) {
      FreePar(i) = par_(freeParIndex(i, p));
      for (int j = 0; j < FreeParCov.GetNrows(); j++) {
        FreeParCov(i, j) = parcov_(freeParIndex(i, p), freeParIndex(j, p));
      }
    }
    TVectorT<double> TrackPar = computePar(FreePar);
    TMatrixTSym<double> TrackCov =
        ErrorMatrixPropagator::propagateError(&TrackHelixVertexFitter::computePar, FreePar, FreeParCov);
    refitParticles.push_back(TrackParticle(TrackPar,
                                           TrackCov,
                                           particles_[p].pdgId(),
                                           particles_[p].mass(),
                                           particles_[p].charge(),
                                           particles_[p].bField()));
  }
  return refitParticles;
}

std::vector<LorentzVectorParticle> TrackHelixVertexFitter::getRefitLorentzVectorParticles() {
  std::vector<LorentzVectorParticle> refitParticles;
  for (unsigned int p = 0; p < particles_.size(); p++) {
    TVectorT<double> FreePar(NFreeTrackPar + NExtraPar + MassOffSet);
    TMatrixTSym<double> FreeParCov(NFreeTrackPar + NExtraPar + MassOffSet);
    for (int i = 0; i < NFreeTrackPar; i++) {
      FreePar(i) = par_(freeParIndex(i, p));
      for (int j = 0; j < NFreeTrackPar; j++) {
        FreeParCov(i, j) = parcov_(freeParIndex(i, p), freeParIndex(j, p));
      }
    }
    FreePar(NFreeTrackPar + MassOffSet) = particles_[p].mass();
    FreePar(NFreeTrackPar + BField0) = particles_[p].bField();
    TVectorT<double> LVPar = computeLorentzVectorPar(FreePar);
    TMatrixTSym<double> LVCov =
        ErrorMatrixPropagator::propagateError(&TrackHelixVertexFitter::computeLorentzVectorPar, FreePar, FreeParCov);
    refitParticles.push_back(
        LorentzVectorParticle(LVPar, LVCov, particles_[p].pdgId(), particles_[p].charge(), particles_[p].bField()));
  }
  return refitParticles;
}

LorentzVectorParticle TrackHelixVertexFitter::getMother(int pdgid) {
  double c(0), b(0);
  TVectorT<double> FreePar(par_.GetNrows() + NExtraPar + particles_.size());
  TMatrixTSym<double> FreeParCov(par_.GetNrows() + NExtraPar + particles_.size());
  for (int i = 0; i < par_.GetNrows(); i++) {
    FreePar(i) = par_(i);
    for (int j = 0; j < par_.GetNrows(); j++) {
      FreeParCov(i, j) = parcov_(i, j);
    }
  }
  for (unsigned int p = 0; p < particles_.size(); p++) {
    b = particles_[p].bField();
    c += particles_[p].charge();
    FreePar(par_.GetNrows() + MassOffSet + p) = particles_[p].mass();
  }
  FreePar(par_.GetNrows() + BField0) = b;
  TVectorT<double> mpar = computeMotherLorentzVectorPar(FreePar);
  TMatrixTSym<double> mcov = ErrorMatrixPropagator::propagateError(
      &TrackHelixVertexFitter::computeMotherLorentzVectorPar, FreePar, FreeParCov);
  return LorentzVectorParticle(mpar, mcov, pdgid, c, b);
}

TVector3 TrackHelixVertexFitter::getVertex() {
  return TVector3(par_(freeParIndex(x0, 0)), par_(freeParIndex(y0, 0)), par_(freeParIndex(z0, 0)));
}

TMatrixTSym<double> TrackHelixVertexFitter::getVertexError() {
  TMatrixTSym<double> c(NFreeVertexPar);
  for (unsigned int i = 0; i < NFreeVertexPar; i++) {
    for (unsigned int j = 0; j < NFreeVertexPar; j++) {
      c(freeParIndex(i, 0), freeParIndex(j, 0)) = parcov_(freeParIndex(i, 0), freeParIndex(j, 0));
    }
  }
  return c;
}

void TrackHelixVertexFitter::computedxydz(const TVectorT<double>& inpar,
                                          int p,
                                          double& kappa,
                                          double& lam,
                                          double& phi,
                                          double& x,
                                          double& y,
                                          double& z,
                                          double& s,
                                          double& dxy,
                                          double& dz) {
  kappa = inpar(freeParIndex(kappa0, p));
  lam = inpar(freeParIndex(lambda0, p));
  phi = inpar(freeParIndex(phi0, p));
  x = inpar(freeParIndex(x0, p));
  y = inpar(freeParIndex(y0, p));
  z = inpar(freeParIndex(z0, p));
  double v = (2.0 * kappa * (x * cos(phi) + y * sin(phi)));
  double arcsinv = 0;
  if (v >= 1.0) {
    arcsinv = TMath::Pi() / 2;
  } else if (v <= -1.0) {
    arcsinv = -TMath::Pi() / 2;
  } else {
    arcsinv = asin(v);
  }
  s = 1.0 / (2.0 * kappa) * arcsinv;
  dxy = y * cos(phi) - x * sin(phi) - (1 / kappa) * sin(kappa * s) * sin(kappa * s);
  dz = z - s * tan(lam);
}

TVectorT<double> TrackHelixVertexFitter::computePar(const TVectorT<double>& inpar) {
  int nparticles = (inpar.GetNrows() - NFreeVertexPar) / (NFreeTrackPar - NFreeVertexPar);
  TVectorT<double> helices(nparticles * TrackParticle::NHelixPar);
  for (int p = 0; p < nparticles; p++) {
    TVectorT<double> TrackPar = computeTrackPar(inpar, p);
    for (int i = 0; i < TrackParticle::NHelixPar; i++) {
      helices(measuredValueIndex(i, p)) = TrackPar(i);
    }
  }
  return helices;
}

TVectorT<double> TrackHelixVertexFitter::computeTrackPar(const TVectorT<double>& inpar, int p) {
  TVectorT<double> helix(TrackParticle::NHelixPar);
  // copy parameters that are 1 to 1
  double kappa, lam, phi, x, y, z, s, dxy, dz;
  TrackHelixVertexFitter::computedxydz(inpar, p, kappa, lam, phi, x, y, z, s, dxy, dz);
  helix(TrackParticle::kappa) = kappa;
  helix(TrackParticle::lambda) = lam;
  helix(TrackParticle::phi) = phi;
  helix(TrackParticle::dxy) = dxy;
  helix(TrackParticle::dz) = dz;
  return helix;
}

TVectorT<double> TrackHelixVertexFitter::computeLorentzVectorPar(const TVectorT<double>& inpar) {
  int np(0), parsize(0);
  parSizeInfo(inpar, np, parsize, true);
  double B = inpar(parsize + BField0);
  double massHypothesis = inpar(parsize + MassOffSet);
  TVectorT<double> LV(LorentzVectorParticle::NLorentzandVertexPar);
  double kappa, lam, phi, x, y, z, s, dxy, dz;
  int p = 0;
  TrackHelixVertexFitter::computedxydz(inpar, p, kappa, lam, phi, x, y, z, s, dxy, dz);
  double phi1 = 2 * s * kappa + phi;
  double bOverK = B * (1.0 / fabs(kappa));
  LV(LorentzVectorParticle::px) = bOverK * cos(phi1);
  LV(LorentzVectorParticle::py) = bOverK * sin(phi1);
  LV(LorentzVectorParticle::pz) = bOverK * tan(lam);
  LV(LorentzVectorParticle::m) = massHypothesis;
  LV(LorentzVectorParticle::vx) = x;
  LV(LorentzVectorParticle::vy) = y;
  LV(LorentzVectorParticle::vz) = z;
  return LV;
}

TVectorT<double> TrackHelixVertexFitter::computeMotherLorentzVectorPar(const TVectorT<double>& inpar) {
  TVectorT<double> mother(LorentzVectorParticle::NLorentzandVertexPar);
  double E(0);
  int np(0), parsize(0);
  parSizeInfo(inpar, np, parsize, true);
  for (int p = 0; p < np; p++) {
    TVectorT<double> particlepar(NFreeTrackPar + NExtraPar + MassOffSet);
    for (int i = 0; i < NFreeTrackPar; i++) {
      particlepar(i) = inpar(freeParIndex(i, p));
    }
    particlepar(NFreeTrackPar + BField0) = inpar(parsize + BField0);
    particlepar(NFreeTrackPar + MassOffSet) = inpar(parsize + MassOffSet + p);
    TVectorT<double> daughter = TrackHelixVertexFitter::computeLorentzVectorPar(particlepar);
    mother(LorentzVectorParticle::px) += daughter(LorentzVectorParticle::px);
    mother(LorentzVectorParticle::py) += daughter(LorentzVectorParticle::py);
    mother(LorentzVectorParticle::pz) += daughter(LorentzVectorParticle::pz);
    mother(LorentzVectorParticle::vx) = daughter(LorentzVectorParticle::vx);
    mother(LorentzVectorParticle::vy) = daughter(LorentzVectorParticle::vy);
    mother(LorentzVectorParticle::vz) = daughter(LorentzVectorParticle::vz);
    E += sqrt((daughter(LorentzVectorParticle::px) * daughter(LorentzVectorParticle::px) +
               daughter(LorentzVectorParticle::py) * daughter(LorentzVectorParticle::py) +
               daughter(LorentzVectorParticle::pz) * daughter(LorentzVectorParticle::pz) +
               daughter(LorentzVectorParticle::m) * daughter(LorentzVectorParticle::m)));
  }
  double P2 = (mother(LorentzVectorParticle::px) * mother(LorentzVectorParticle::px) +
               mother(LorentzVectorParticle::py) * mother(LorentzVectorParticle::py) +
               mother(LorentzVectorParticle::pz) * mother(LorentzVectorParticle::pz));
  mother(LorentzVectorParticle::m) = (E * E - P2) / sqrt(fabs(E * E - P2));
  return mother;
}

TString TrackHelixVertexFitter::freeParName(int Par) {
  int p(0);
  if (Par == x0)
    return "x0";
  if (Par == y0)
    return "y0";
  if (Par == z0)
    return "z0";
  for (p = 0; p < nParticles_; p++) {
    if ((Par - NFreeVertexPar) < (p + 1) * (NFreeTrackPar - NFreeVertexPar))
      break;
  }
  TString n;
  int index = Par - p * (NFreeTrackPar - NFreeVertexPar);
  if (index == kappa0)
    n = "kappa0";
  if (index == lambda0)
    n = "lambda0";
  if (index == phi0)
    n = "phi0";
  n += "_particle";
  n += p;
  return n;
}

void TrackHelixVertexFitter::parSizeInfo(const TVectorT<double>& inpar, int& np, int& parsize, bool hasextras) {
  if (hasextras)
    np = (inpar.GetNrows() - NFreeVertexPar - NExtraPar) / (NFreeTrackPar + MassOffSet - NFreeVertexPar);
  else
    np = (inpar.GetNrows() - NFreeVertexPar) / (NFreeTrackPar - NFreeVertexPar);
  parsize = np * (NFreeTrackPar - NFreeVertexPar) + NFreeVertexPar;
}