File indexing completed on 2024-04-06 12:31:33
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
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
0046 }
0047
0048
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
0057 virtual SubDetector subDetector() const { return GeomDetEnumerators::DT; }
0058 };
0059
0060
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
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 }