Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:31:37

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 }