KinematicFitPrimitives/test/KinematicState_t.cpp

0001 #include "RecoVertex/KinematicFitPrimitives/interface/TrackKinematicStatePropagator.h"
0002 #include "RecoVertex/KinematicFitPrimitives/interface/KinematicState.h"
0003 #include "RecoVertex/KinematicFitPrimitives/interface/Matrices.h"
0004
0005 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0006 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateAccessor.h"
0007 #include "DataFormats/GeometrySurface/interface/Surface.h"
0008 #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
0009 #include "MagneticField/Engine/interface/MagneticField.h"
0010
0011 #include "TrackingTools/TrajectoryState/interface/BasicSingleTrajectoryState.h"
0012 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0013 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCurvilinearToCartesian.h"
0014 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCartesianToCurvilinear.h"
0015 #include "TrackingTools/TrajectoryState/interface/PerigeeConversions.h"
0016
0017 #include <iostream>
0018
0019 class ConstMagneticField final : public MagneticField {
0020 public:
0021   ConstMagneticField() { setNominalValue(); }
0022
0023   GlobalVector inTesla(const GlobalPoint&) const override { return GlobalVector(0, 0, 4); }
0024 };
0025
0026 int main() {
0027   using namespace std;
0028
0029   MagneticField* field = new ConstMagneticField;
0030   GlobalPoint gp(0, 0, 0);
0031   GlobalPoint gp2(1, 1, 1);
0032
0033   GlobalPoint stars(-1.e12, -1.e12, -1);
0034
0035   GlobalVector gv(1, 1, 1);
0036   GlobalTrajectoryParameters gtp(gp, gv, 1, field);
0037   double v[15] = {0.01, -0.01, 0., 0., 0., 0.01, 0., 0., 0., 0.01, 0., 0., 1., 0., 1.};
0038   AlgebraicSymMatrix55 gerr(v, 15);
0039   BoundPlane* plane = new BoundPlane(gp, Surface::RotationType());
0040
0041   TrajectoryStateOnSurface ts(gtp, gerr, *plane);
0042
0043   cout << "ts.globalPosition() " << ts.globalPosition() << endl;
0044   cout << "ts.globalMomentum() " << ts.globalMomentum() << endl;
0045   cout << "ts.localMomentum()  " << ts.localMomentum() << endl;
0046   cout << "ts.transverseCurvature()  " << ts.transverseCurvature() << endl;
0047   cout << "ts inversePtErr " << TrajectoryStateAccessor(*ts.freeState()).inversePtError() << std::endl;
0048   cout << "ts curv err\n" << ts.curvilinearError().matrix() << std::endl;
0049   cout << "ts cart err\n" << ts.cartesianError().matrix() << std::endl;
0050
0051   const FreeTrajectoryState* fts = ts.freeTrajectoryState();
0052   assert(fts);
0053
0054   KinematicState ks(*fts, 0.151, 0.01);
0055
0056   std::cout << "ks Par" << ks.kinematicParameters().vector() << std::endl;
0057   std::cout << "ks Err" << ks.kinematicParametersError().matrix() << std::endl;
0058   std::cout << "ks charge " << ks.particleCharge() << std::endl;
0059   cout << "ks.globalPosition() " << ks.globalPosition() << endl;
0060   cout << "ks globalMomentum() " << ks.globalMomentum() << endl;
0061   cout << "ks curv err\n" << ks.freeTrajectoryState().curvilinearError().matrix() << std::endl;
0062   cout << "ks cart err\n" << ks.freeTrajectoryState().cartesianError().matrix() << std::endl;
0063
0064   auto c = ks.kinematicParametersError().matrix();
0065   c(0, 6) = -0.0001;
0066   KinematicParametersError kpe(c);
0067   KinematicState ks2(ks.kinematicParameters(), kpe, ks.particleCharge(), ks.magneticField());
0068
0069   std::cout << "ks2 Par" << ks2.kinematicParameters().vector() << std::endl;
0070   std::cout << "ks2 Err" << ks2.kinematicParametersError().matrix() << std::endl;
0071   std::cout << "ks2 charge " << ks2.particleCharge() << std::endl;
0072   cout << "ks2.globalPosition() " << ks2.globalPosition() << endl;
0073   cout << "ks2 globalMomentum() " << ks2.globalMomentum() << endl;
0074   cout << "ks2 curv err\n" << ks2.freeTrajectoryState().curvilinearError().matrix() << std::endl;
0075   cout << "ks2 cart err\n" << ks2.freeTrajectoryState().cartesianError().matrix() << std::endl;
0076
0077   TrackKinematicStatePropagator p;
0078   {
0079     auto ok = p.willPropagateToTheTransversePCA(ks, stars);
0080     std::cout << "\npropagate ks to stars " << (ok ? "ok\n" : "nope\n") << std::endl;
0081   }
0082
0083   auto ok = p.willPropagateToTheTransversePCA(ks, gp2);
0084   std::cout << "\npropagate ks " << (ok ? "ok\n" : "nope\n") << std::endl;
0085   auto kst = p.propagateToTheTransversePCA(ks, gp2);
0086
0087   std::cout << "kst Par" << kst.kinematicParameters().vector() << std::endl;
0088   std::cout << "kst Err" << kst.kinematicParametersError().matrix() << std::endl;
0089   cout << "kst.globalPosition() " << kst.globalPosition() << endl;
0090   cout << "kst globalMomentum() " << kst.globalMomentum() << endl;
0091   cout << "kst curv err\n" << kst.freeTrajectoryState().curvilinearError().matrix() << std::endl;
0092   cout << "kst cart err\n" << kst.freeTrajectoryState().cartesianError().matrix() << std::endl;
0093
0094   auto ok2 = p.willPropagateToTheTransversePCA(ks2, gp2);
0095   std::cout << "\npropagate ks2 " << (ok2 ? "ok\n" : "nope\n") << std::endl;
0096   auto kst2 = p.propagateToTheTransversePCA(ks2, gp2);
0097
0098   std::cout << "kst2 Par" << kst2.kinematicParameters().vector() << std::endl;
0099   std::cout << "kst2 Err" << kst2.kinematicParametersError().matrix() << std::endl;
0100   cout << "kst2.globalPosition() " << kst2.globalPosition() << endl;
0101   cout << "kst2 globalMomentum() " << kst2.globalMomentum() << endl;
0102   cout << "kst2 curv err\n" << kst2.freeTrajectoryState().curvilinearError().matrix() << std::endl;
0103   cout << "kst2 cart err\n" << kst2.freeTrajectoryState().cartesianError().matrix() << std::endl;
0104
0105   return 0;
0106 }