ReferenceTrajectories/src/TwoBodyDecayTrajectoryState.cc

0001
0002 #include "Alignment/ReferenceTrajectories/interface/TwoBodyDecayTrajectoryState.h"
0003 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0004
0005 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayModel.h"
0006 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayDerivatives.h"
0007
0008 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0009
0010 #include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h"
0011 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCurvilinearToLocal.h"
0012 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCartesianToCurvilinear.h"
0013
0014 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0015
0016 #include "DataFormats/CLHEP/interface/Migration.h"
0017
0018 using namespace std;
0019
0020 TwoBodyDecayTrajectoryState::TwoBodyDecayTrajectoryState(const TsosContainer& tsos,
0021                                                          const TwoBodyDecay& tbd,
0022                                                          double particleMass,
0023                                                          const MagneticField* magField,
0024                                                          bool propagateErrors)
0025     : theValidityFlag(false),
0026       theParticleMass(particleMass),
0027       theParameters(tbd.decayParameters()),
0028       theDerivatives(AlgebraicMatrix(nLocalParam, nDecayParam), AlgebraicMatrix(nLocalParam, nDecayParam)),
0029       theOriginalTsos(tsos),
0030       thePrimaryMass(tbd.primaryMass()),
0031       thePrimaryWidth(tbd.primaryWidth()) {
0032   construct(magField, propagateErrors);
0033 }
0034
0035 void TwoBodyDecayTrajectoryState::rescaleError(double scale) {
0036   theOriginalTsos.first.rescaleError(scale);
0037   theOriginalTsos.second.rescaleError(scale);
0038   theRefittedTsos.first.rescaleError(scale);
0039   theRefittedTsos.second.rescaleError(scale);
0040 }
0041
0042 void TwoBodyDecayTrajectoryState::construct(const MagneticField* magField, bool propagateErrors) {
0043   // construct global trajectory parameters at the starting point
0044   TwoBodyDecayModel tbdDecayModel(theParameters[TwoBodyDecayParameters::mass], theParticleMass);
0045   pair<AlgebraicVector, AlgebraicVector> secondaryMomenta = tbdDecayModel.cartesianSecondaryMomenta(theParameters);
0046
0047   GlobalPoint vtx(theParameters[TwoBodyDecayParameters::x],
0048                   theParameters[TwoBodyDecayParameters::y],
0049                   theParameters[TwoBodyDecayParameters::z]);
0050
0051   GlobalVector p1(secondaryMomenta.first[0], secondaryMomenta.first[1], secondaryMomenta.first[2]);
0052
0053   GlobalVector p2(secondaryMomenta.second[0], secondaryMomenta.second[1], secondaryMomenta.second[2]);
0054
0055   GlobalTrajectoryParameters gtp1(vtx, p1, theOriginalTsos.first.charge(), magField);
0056   FreeTrajectoryState fts1(gtp1);
0057
0058   GlobalTrajectoryParameters gtp2(vtx, p2, theOriginalTsos.second.charge(), magField);
0059   FreeTrajectoryState fts2(gtp2);
0060
0061   // contruct derivatives at the starting point
0062   TwoBodyDecayDerivatives tbdDerivatives(theParameters[TwoBodyDecayParameters::mass], theParticleMass);
0063   pair<AlgebraicMatrix, AlgebraicMatrix> derivatives = tbdDerivatives.derivatives(theParameters);
0064
0065   AlgebraicMatrix deriv1(6, 9, 0);
0066   deriv1.sub(1, 1, AlgebraicMatrix(3, 3, 1));
0067   deriv1.sub(4, 4, derivatives.first);
0068
0069   AlgebraicMatrix deriv2(6, 9, 0);
0070   deriv2.sub(1, 1, AlgebraicMatrix(3, 3, 1));
0071   deriv2.sub(4, 4, derivatives.second);
0072
0073   // compute errors of initial states
0074   if (propagateErrors) {
0075     setError(fts1, deriv1);
0076     setError(fts2, deriv2);
0077   }
0078
0079   // propgate states and derivatives from the starting points to the end points
0080   bool valid1 = propagateSingleState(
0081       fts1, gtp1, deriv1, theOriginalTsos.first.surface(), magField, theRefittedTsos.first, theDerivatives.first);
0082
0083   bool valid2 = propagateSingleState(
0084       fts2, gtp2, deriv2, theOriginalTsos.second.surface(), magField, theRefittedTsos.second, theDerivatives.second);
0085
0086   theValidityFlag = valid1 && valid2;
0087
0088   return;
0089 }
0090
0091 bool TwoBodyDecayTrajectoryState::propagateSingleState(const FreeTrajectoryState& fts,
0092                                                        const GlobalTrajectoryParameters& gtp,
0093                                                        const AlgebraicMatrix& startDeriv,
0094                                                        const Surface& surface,
0095                                                        const MagneticField* magField,
0096                                                        TrajectoryStateOnSurface& tsos,
0097                                                        AlgebraicMatrix& endDeriv) const {
0098   AnalyticalPropagator propagator(magField);
0099
0100   // propagate state
0101   pair<TrajectoryStateOnSurface, double> tsosWithPath = propagator.propagateWithPath(fts, surface);
0102
0103   // check if propagation was successful
0104   if (!tsosWithPath.first.isValid())
0105     return false;
0106
0107   // jacobian for transformation from cartesian to curvilinear frame at the starting point
0108   JacobianCartesianToCurvilinear cartToCurv(gtp);
0109   const AlgebraicMatrix56& matCartToCurv = cartToCurv.jacobian();
0110
0111   // jacobian in curvilinear frame for propagation from the starting point to the end point
0112   AnalyticalCurvilinearJacobian curvJac(
0113       gtp, tsosWithPath.first.globalPosition(), tsosWithPath.first.globalMomentum(), tsosWithPath.second);
0114   const AlgebraicMatrix55& matCurvJac = curvJac.jacobian();
0115
0116   // jacobian for transformation from curvilinear to local frame at the end point
0117   JacobianCurvilinearToLocal curvToLoc(surface, tsosWithPath.first.localParameters(), *magField);
0118   const AlgebraicMatrix55& matCurvToLoc = curvToLoc.jacobian();
0119
0120   AlgebraicMatrix56 tmpDeriv = matCurvToLoc * matCurvJac * matCartToCurv;
0121   AlgebraicMatrix hepMatDeriv(asHepMatrix(tmpDeriv));
0122   //AlgebraicMatrix hepMatDeriv = asHepMatrix< 5, 6 >( tmpDeriv );
0123
0124   // replace original state with new state
0125   tsos = tsosWithPath.first;
0126
0127   // propagate derivative matrix
0128   endDeriv = hepMatDeriv * startDeriv;
0129
0130   return true;
0131 }
0132
0133 void TwoBodyDecayTrajectoryState::setError(FreeTrajectoryState& fts, AlgebraicMatrix& derivative) const {
0134   AlgebraicSymMatrix ftsCartesianError(theParameters.covariance().similarity(derivative));
0135   fts.setCartesianError(asSMatrix<6>(ftsCartesianError));
0136 }