GsfTools/test/traceSymMult_t.cpp

0001 #include "TrackingTools/GsfTools/interface/KullbackLeiblerDistance.h"
0002 #include "TrackingTools/GsfTools/interface/DistanceBetweenComponents.h"
0003
0004 #include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCartesian.h"
0005 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCartesianToLocal.h"
0006 #include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCurvilinear.h"
0007
0008 #include "DataFormats/GeometrySurface/interface/Surface.h"
0009 #include "TrackingTools/TrajectoryParametrization/interface/LocalTrajectoryParameters.h"
0010
0011 #include "DataFormats/GeometrySurface/interface/Plane.h"
0012
0013 #include "TrackingTools/GsfTools/interface/GsfMatrixTools.h"
0014
0015 #include "FWCore/Utilities/interface/HRRealTime.h"
0016 #include <iostream>
0017
0018 void st() {}
0019 void en() {}
0020
0021 typedef DistanceBetweenComponents<5> Distance;
0022 typedef KullbackLeiblerDistance<5> KDistance;
0023 typedef SingleGaussianState<5> __attribute__((aligned(16))) GS;
0024 typedef GS::Vector __attribute__((aligned(16))) Vector;
0025 typedef GS::Matrix __attribute__((aligned(16))) Matrix;
0026 typedef ROOT::Math::SMatrix<double, 6, 6, ROOT::Math::MatRepSym<double, 6> > __attribute__((aligned(16))) Matrix6;
0027 typedef ROOT::Math::SMatrix<double, 2, 2, ROOT::Math::MatRepSym<double, 2> > __attribute__((aligned(16))) Matrix2;
0028
0029 Distance const& distance() {
0030   static Distance* d = new KDistance;
0031   return *d;
0032 }
0033
0034 Matrix buildCovariance(float y) {
0035   // build a resonable covariance matrix as JIJ
0036
0037   Basic3DVector<float> axis(0.5, 1., 1);
0038
0039   Surface::RotationType rot(axis, 0.5 * M_PI);
0040
0041   Surface::PositionType pos(0., 0., 0.);
0042
0043   Plane plane(pos, rot);
0044   LocalTrajectoryParameters tp(1., 1., y, 0., 0., 1.);
0045
0046   JacobianLocalToCartesian jl2c(plane, tp);
0047   return ROOT::Math::SimilarityT(jl2c.jacobian(), Matrix6(ROOT::Math::SMatrixIdentity()));
0048   // return  ROOT::Math::Transpose(jl2c.jacobian())* jl2c.jacobian();
0049 }
0050
0051 int main(int args, char** argv) {
0052   std::cout << "size of 2x2 Sym Matrix " << sizeof(Matrix2) << std::endl;
0053   std::cout << "size of 5x5 Sym Matrix " << sizeof(Matrix) << std::endl;
0054   std::cout << "size of 6x6 Sym Matrix " << sizeof(Matrix6) << std::endl;
0055
0056   // Distance const & d = distance();
0057
0058   Matrix cov1 = buildCovariance(2.);
0059   Matrix cov2 = buildCovariance(1.);
0060
0061   double one = GsfMatrixTools::trace(cov1, cov2);
0062
0063   double two = GsfMatrixTools::trace<double, 5>(cov1 * cov2);
0064
0065   if (fabs(one - two) > 1.e-12)
0066     std::cout << "vincenzo was wrong!" << std::endl;
0067   std::cout << one << " " << two << " " << one - two << std::endl;
0068
0069   if (args == 1) {
0070     st();
0071     one = GsfMatrixTools::trace(cov2, cov1);
0072     en();
0073   } else {
0074     st();
0075     two = GsfMatrixTools::trace<double, 5>(cov2 * cov1);
0076     en();
0077   }
0078
0079   if (fabs(one - two) > 1.e-15)
0080     std::cout << "vincenzo was wrong!" << std::endl;
0081
0082   AlgebraicVector5 v(1., 0.5, -0.5, 0.3, 0.7);
0083   AlgebraicSymMatrix55 c1;
0084   AlgebraicSymMatrix55 c2;
0085   AlgebraicSymMatrix11 s = AlgebraicMatrixID();  //stupid trick to make CLHEP work decently
0086   c1 = ROOT::Math::Similarity(AlgebraicMatrix51(v.Array(), 5), s);
0087   std::cout << c1 << std::endl;
0088   ROOT::Math::AssignSym::Evaluate(c2, ROOT::Math::TensorProd(v, v));
0089   std::cout << c2 << std::endl;
0090
0091   return 0;
0092 }