KalmanUpdators/test/KFUpdator_t.cpp

0001 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0002 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0003
0004 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0005 #include "DataFormats/GeometrySurface/interface/Surface.h"
0006 #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
0007 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
0008
0009 #include "MagneticField/Engine/interface/MagneticField.h"
0010
0011 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0012 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0013 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
0014 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0015 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
0016 #include <iostream>
0017
0018 class ConstMagneticField final : public MagneticField {
0019 public:
0020   ConstMagneticField() { setNominalValue(); }
0021   GlobalVector inTesla(const GlobalPoint&) const override { return GlobalVector(0, 0, 4); }
0022 };
0023
0024 #include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCartesian.h"
0025 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCartesianToLocal.h"
0026 #include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCurvilinear.h"
0027
0028 typedef ROOT::Math::SMatrix<double, 5, 5, ROOT::Math::MatRepSym<double, 5> > Matrix5;
0029 typedef ROOT::Math::SMatrix<double, 6, 6, ROOT::Math::MatRepSym<double, 6> > Matrix6;
0030
0031 Matrix5 buildCovariance(float y) {
0032   // build a resonable covariance matrix as JIJ
0033
0034   Basic3DVector<float> axis(0.5, 1., 1);
0035
0036   Surface::RotationType rot(axis, 0.5 * M_PI);
0037
0038   Surface::PositionType pos(0., 0., 0.);
0039
0040   Plane plane(pos, rot);
0041   LocalTrajectoryParameters tp(1., 1., y, 0., 0., 1.);
0042
0043   JacobianLocalToCartesian jl2c(plane, tp);
0044   return ROOT::Math::SimilarityT(jl2c.jacobian(), Matrix6(ROOT::Math::SMatrixIdentity()));
0045   // return  ROOT::Math::Transpose(jl2c.jacobian())* jl2c.jacobian();
0046 }
0047
0048 // A fake Det class
0049
0050 class MyDet : public TrackerGeomDet {
0051 public:
0052   MyDet(BoundPlane* bp, DetId id) : TrackerGeomDet(bp) { setDetId(id); }
0053
0054   virtual std::vector<const GeomDet*> components() const { return std::vector<const GeomDet*>(); }
0055
0056   /// Which subdetector
0057   virtual SubDetector subDetector() const { return GeomDetEnumerators::DT; }
0058 };
0059
0060 // a 2d hit
0061
0062 #include "TrackingTools/TransientTrackingRecHit/interface/GenericTransientTrackingRecHit.h"
0063 #include "TrackingTools/TransientTrackingRecHit/interface/HelpertRecHit2DLocalPos.h"
0064
0065 class My2DHit : public GenericTransientTrackingRecHit {
0066 public:
0067   My2DHit(const GeomDet& geom, TrackingRecHit* rh) : GenericTransientTrackingRecHit(geom, rh) {}
0068
0069   virtual void getKfComponents(KfComponentsHolder& holder) const {
0070     HelpertRecHit2DLocalPos::getKfComponents(holder, *hit(), *det());
0071   }
0072 };
0073
0074 void print(TrajectoryStateOnSurface const& ts) {
0075   using namespace std;
0076
0077   cout << "transverseCurvature  " << ts.transverseCurvature() << endl;
0078   cout << "globalMomentum       " << ts.globalMomentum() << endl;
0079   cout << "localMomentum        " << ts.localMomentum() << endl;
0080   cout << "localPosition        " << ts.localPosition() << endl;
0081   cout << "localError           " << ts.localError().matrix() << endl;
0082   cout << endl;
0083 }
0084
0085 #include "FWCore/Utilities/interface/HRRealTime.h"
0086 #include <iostream>
0087 #include <vector>
0088
0089 bool isAligned(const void* data, long alignment) {
0090   // check that the alignment is a power of two
0091   assert((alignment & (alignment - 1)) == 0);
0092   return ((long)data & (alignment - 1)) == 0;
0093 }
0094
0095 void st() {}
0096 void en() {}
0097
0098 struct KFUTest {
0099   TrajectoryStateUpdator const& tsu;
0100
0101   KFUTest(TrajectoryStateUpdator const* itsu) : tsu(*itsu) {}
0102
0103   void print(const TrajectoryStateOnSurface& tsos, const TrackingRecHit& hit) const {
0104     TrajectoryStateOnSurface tsn = tsu.update(tsos, hit);
0105     ::print(tsn);
0106   }
0107
0108   void time(const TrajectoryStateOnSurface& tsos, const TrackingRecHit& hit) const {
0109     edm::HRTimeType s = edm::hrRealTime();
0110     st();
0111     TrajectoryStateOnSurface tsn = tsu.update(tsos, hit);
0112     en();
0113     edm::HRTimeType e = edm::hrRealTime();
0114     std::cout << e - s << std::endl;
0115   }
0116 };
0117
0118 struct Chi2Test {
0119   MeasurementEstimator const& chi2;
0120   mutable std::pair<bool, double> res;
0121   Chi2Test(MeasurementEstimator const* ichi2) : chi2(*ichi2) {}
0122
0123   void print(const TrajectoryStateOnSurface& tsos, const TrackingRecHit& hit) const {
0124     res = chi2.estimate(tsos, hit);
0125     std::cout << "chi2 " << res.second << std::endl;
0126   }
0127
0128   void time(const TrajectoryStateOnSurface& tsos, const TrackingRecHit& hit) const {
0129     edm::HRTimeType s = edm::hrRealTime();
0130     st();
0131     res = chi2.estimate(tsos, hit);
0132     en();
0133     edm::HRTimeType e = edm::hrRealTime();
0134     std::cout << e - s << std::endl;
0135   }
0136 };
0137
0138 int main() {
0139   MagneticField* field = new ConstMagneticField;
0140   GlobalPoint gp(0, 0, 0);
0141   GlobalVector gv(1, 1, 1);
0142   GlobalTrajectoryParameters gtp(gp, gv, 1, field);
0143   CurvilinearTrajectoryError gerr(buildCovariance(1.));
0144   BoundPlane* plane = new BoundPlane(gp, Surface::RotationType());
0145   GeomDet* det = new MyDet(plane, 41);
0146   GeomDet* gdet = new MyDet(plane, 40);
0147
0148   TrajectoryStateOnSurface ts(gtp, gerr, *plane);
0149   print(ts);
0150
0151   LocalPoint lp(0, 0, 0);
0152   LocalVector lv(1, 1, 1);
0153   LocalTrajectoryParameters ltp(lp, lv, 1);
0154   LocalTrajectoryError ler(0.1, 0.1, 0.01, 0.05, 0.1);
0155   TrajectoryStateOnSurface ts2(ltp, ler, *plane, field);
0156   print(ts2);
0157
0158   LocalPoint m(0.1, 0.1, 0);
0159   LocalError e(0.2, -0.05, 0.1);
0160
0161   OmniClusterRef cref;
0162   SiPixelRecHit::ClusterRef pref;
0163   SiStripRecHit2D dummy;
0164   SiPixelRecHit hitpx(m, e, 1., *det, pref);
0165   SiStripRecHit1D hit1d(m, e, *det, cref);
0166   SiStripRecHit2D hit2d(m, e, *det, cref);
0167   ProjectedSiStripRecHit2D hitpj(m, e, *gdet, hit2d);
0168   TrackingRecHit* hit = new SiStripMatchedRecHit2D(m, e, *det, &dummy, &dummy);
0169   TransientTrackingRecHit* thit = new My2DHit(*det, hit);
0170
0171   KFUTest kt(new KFUpdator());
0172   Chi2Test chi2(new Chi2MeasurementEstimator(10.));
0173
0174   std::cout << "\n** KFU ** \n" << std::endl;
0175
0176   kt.print(ts, *thit);
0177   kt.print(ts2, *thit);
0178
0179   kt.time(ts, *thit);
0180   kt.time(ts2, *thit);
0181
0182   kt.print(ts, hit2d);
0183   kt.print(ts2, hit2d);
0184
0185   kt.print(ts, hitpx);
0186   kt.print(ts2, hitpx);
0187
0188   kt.print(ts, hitpj);
0189   kt.print(ts2, hitpj);
0190
0191   kt.print(ts, hit1d);
0192   kt.print(ts2, hit1d);
0193
0194   kt.time(ts, *thit);
0195   kt.time(ts2, *thit);
0196
0197   std::cout << "\n** Chi2 ** \n" << std::endl;
0198
0199   chi2.print(ts, *thit);
0200   chi2.print(ts2, *thit);
0201
0202   chi2.print(ts, hit2d);
0203   chi2.print(ts2, hit2d);
0204
0205   chi2.print(ts, hitpx);
0206   chi2.print(ts2, hitpx);
0207
0208   chi2.print(ts, hitpj);
0209   chi2.print(ts2, hitpj);
0210
0211   chi2.print(ts, hit1d);
0212   chi2.print(ts2, hit1d);
0213
0214   chi2.time(ts, *thit);
0215   chi2.time(ts2, *thit);
0216
0217   return 0;
0218 }