TwoBodyDecayTrajectoryState.cc

Go to the documentation of this file.

#include "Alignment/ReferenceTrajectories/interface/TwoBodyDecayTrajectoryState.h"
#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"

#include "Alignment/TwoBodyDecay/interface/TwoBodyDecayModel.h"
#include "Alignment/TwoBodyDecay/interface/TwoBodyDecayDerivatives.h"

#include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"

#include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h"
#include "TrackingTools/AnalyticalJacobians/interface/JacobianCurvilinearToLocal.h"
#include "TrackingTools/AnalyticalJacobians/interface/JacobianCartesianToCurvilinear.h"

#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"

#include "DataFormats/CLHEP/interface/Migration.h"

using namespace std;

TwoBodyDecayTrajectoryState::TwoBodyDecayTrajectoryState(const TsosContainer& tsos,
                                                         const TwoBodyDecay& tbd,
                                                         double particleMass,
                                                         const MagneticField* magField,
                                                         bool propagateErrors)
    : theValidityFlag(false),
      theParticleMass(particleMass),
      theParameters(tbd.decayParameters()),
      theDerivatives(AlgebraicMatrix(nLocalParam, nDecayParam), AlgebraicMatrix(nLocalParam, nDecayParam)),
      theOriginalTsos(tsos),
      thePrimaryMass(tbd.primaryMass()),
      thePrimaryWidth(tbd.primaryWidth()) {
  construct(magField, propagateErrors);
}

void TwoBodyDecayTrajectoryState::rescaleError(double scale) {
  theOriginalTsos.first.rescaleError(scale);
  theOriginalTsos.second.rescaleError(scale);
  theRefittedTsos.first.rescaleError(scale);
  theRefittedTsos.second.rescaleError(scale);
}

void TwoBodyDecayTrajectoryState::construct(const MagneticField* magField, bool propagateErrors) {
  // construct global trajectory parameters at the starting point
  TwoBodyDecayModel tbdDecayModel(theParameters[TwoBodyDecayParameters::mass], theParticleMass);
  pair<AlgebraicVector, AlgebraicVector> secondaryMomenta = tbdDecayModel.cartesianSecondaryMomenta(theParameters);

  GlobalPoint vtx(theParameters[TwoBodyDecayParameters::x],
                  theParameters[TwoBodyDecayParameters::y],
                  theParameters[TwoBodyDecayParameters::z]);

  GlobalVector p1(secondaryMomenta.first[0], secondaryMomenta.first[1], secondaryMomenta.first[2]);

  GlobalVector p2(secondaryMomenta.second[0], secondaryMomenta.second[1], secondaryMomenta.second[2]);

  GlobalTrajectoryParameters gtp1(vtx, p1, theOriginalTsos.first.charge(), magField);
  FreeTrajectoryState fts1(gtp1);

  GlobalTrajectoryParameters gtp2(vtx, p2, theOriginalTsos.second.charge(), magField);
  FreeTrajectoryState fts2(gtp2);

  // contruct derivatives at the starting point
  TwoBodyDecayDerivatives tbdDerivatives(theParameters[TwoBodyDecayParameters::mass], theParticleMass);
  pair<AlgebraicMatrix, AlgebraicMatrix> derivatives = tbdDerivatives.derivatives(theParameters);

  AlgebraicMatrix deriv1(6, 9, 0);
  deriv1.sub(1, 1, AlgebraicMatrix(3, 3, 1));
  deriv1.sub(4, 4, derivatives.first);

  AlgebraicMatrix deriv2(6, 9, 0);
  deriv2.sub(1, 1, AlgebraicMatrix(3, 3, 1));
  deriv2.sub(4, 4, derivatives.second);

  // compute errors of initial states
  if (propagateErrors) {
    setError(fts1, deriv1);
    setError(fts2, deriv2);
  }

  // propgate states and derivatives from the starting points to the end points
  bool valid1 = propagateSingleState(
      fts1, gtp1, deriv1, theOriginalTsos.first.surface(), magField, theRefittedTsos.first, theDerivatives.first);

  bool valid2 = propagateSingleState(
      fts2, gtp2, deriv2, theOriginalTsos.second.surface(), magField, theRefittedTsos.second, theDerivatives.second);

  theValidityFlag = valid1 && valid2;

  return;
}

bool TwoBodyDecayTrajectoryState::propagateSingleState(const FreeTrajectoryState& fts,
                                                       const GlobalTrajectoryParameters& gtp,
                                                       const AlgebraicMatrix& startDeriv,
                                                       const Surface& surface,
                                                       const MagneticField* magField,
                                                       TrajectoryStateOnSurface& tsos,
                                                       AlgebraicMatrix& endDeriv) const {
  AnalyticalPropagator propagator(magField);

  // propagate state
  pair<TrajectoryStateOnSurface, double> tsosWithPath = propagator.propagateWithPath(fts, surface);

  // check if propagation was successful
  if (!tsosWithPath.first.isValid())
    return false;

  // jacobian for transformation from cartesian to curvilinear frame at the starting point
  JacobianCartesianToCurvilinear cartToCurv(gtp);
  const AlgebraicMatrix56& matCartToCurv = cartToCurv.jacobian();

  // jacobian in curvilinear frame for propagation from the starting point to the end point
  AnalyticalCurvilinearJacobian curvJac(
      gtp, tsosWithPath.first.globalPosition(), tsosWithPath.first.globalMomentum(), tsosWithPath.second);
  const AlgebraicMatrix55& matCurvJac = curvJac.jacobian();

  // jacobian for transformation from curvilinear to local frame at the end point
  JacobianCurvilinearToLocal curvToLoc(surface, tsosWithPath.first.localParameters(), *magField);
  const AlgebraicMatrix55& matCurvToLoc = curvToLoc.jacobian();

  AlgebraicMatrix56 tmpDeriv = matCurvToLoc * matCurvJac * matCartToCurv;
  AlgebraicMatrix hepMatDeriv(asHepMatrix(tmpDeriv));
  //AlgebraicMatrix hepMatDeriv = asHepMatrix< 5, 6 >( tmpDeriv );

  // replace original state with new state
  tsos = tsosWithPath.first;

  // propagate derivative matrix
  endDeriv = hepMatDeriv * startDeriv;

  return true;
}

void TwoBodyDecayTrajectoryState::setError(FreeTrajectoryState& fts, AlgebraicMatrix& derivative) const {
  AlgebraicSymMatrix ftsCartesianError(theParameters.covariance().similarity(derivative));
  fts.setCartesianError(asSMatrix<6>(ftsCartesianError));
}