10#include <gtest/gtest.h>
12#include <analysis/VertexFitting/TreeFitter/HelixUtils.h>
13#include <framework/dataobjects/Helix.h>
14#include <framework/geometry/BFieldManager.h>
15#include <framework/gearbox/Const.h>
20 class TreeFitterHelixJacobianTest :
public ::testing::Test {
25 TEST_F(TreeFitterHelixJacobianTest, Helix)
27 const double eps = 1
E-6;
29 const double px = 0.123214;
30 const double py = 1.543;
31 const double pz = -2.876;
32 const double x = 3.765346;
33 const double y = -2.56756;
34 const double z = 1.678;
39 const double aq =
charge / alpha;
43 const double pt = std::sqrt(px * px + py * py);
44 const double omega = aq / pt;
46 const double phi = atan(py / px);
47 const double cosPhi = std::cos(phi);
48 const double sinPhi = std::sin(phi);
50 const double para = -x * cosPhi - y * sinPhi;
51 const double ortho = -y * cosPhi + x * sinPhi;
52 const double R2 = para * para + ortho * ortho;
53 const double A = 2 * ortho + omega * R2;
54 const double U = std::sqrt(1 + omega * A);
55 const double d0 = A / (1 + U);
56 const double l = 1 / omega * atan((omega * para) / (1 + omega * ortho));
57 const double phi0 = phi + atan((omega * para) / (1 + omega * ortho));
58 const double tanLambda = pz / pt;
59 const double z0 = z + tanLambda * l;
61 std::vector<double> h = {d0, phi0, omega, z0, tanLambda};
62 std::vector<double> h_framework = {helix.getD0(), helix.getPhi0(), helix.getOmega(), helix.getZ0(), helix.getTanLambda()};
63 for (
int row = 0; row < 5; ++row) {
64 double res = h[row] - h_framework[row];
65 EXPECT_TRUE(res < eps) <<
"row " << row <<
" num - ana " << res <<
" framework " << h_framework[row] <<
" mine " << h[row];
70 TEST_F(TreeFitterHelixJacobianTest, Parameters)
72 const double delta = 1e-6;
73 const double eps = 1e-5;
75 Eigen::Matrix<double, 5, 6> jacobian_numerical = Eigen::Matrix<double, 5, 6>::Zero(5, 6);
76 Eigen::Matrix<double, 5, 6> jacobian_analytical = Eigen::Matrix<double, 5, 6>::Zero(5, 6);
78 ROOT::Math::XYZVector postmp;
79 ROOT::Math::XYZVector momtmp;
82 const double px = 0.523214;
83 const double py = -1.543;
84 const double pz = -2.876;
85 const double x = -3.765346;
86 const double y = -2.56756;
87 const double z = 5.678;
89 const Eigen::Matrix<double, 1, 6> positionAndMom_ = (Eigen::Matrix<double, 1, 6>() << x, y, z, px, py, pz).finished();
95 for (
int jin = 0; jin < 6; ++jin) {
96 postmp.SetCoordinates(positionAndMom_(0), positionAndMom_(1), positionAndMom_(2));
97 momtmp.SetCoordinates(positionAndMom_(3), positionAndMom_(4), positionAndMom_(5));
98 if (jin == 0) postmp.SetX(postmp.X() + delta);
99 if (jin == 1) postmp.SetY(postmp.Y() + delta);
100 if (jin == 2) postmp.SetZ(postmp.Z() + delta);
101 if (jin == 3) momtmp.SetX(momtmp.X() + delta);
102 if (jin == 4) momtmp.SetY(momtmp.Y() + delta);
103 if (jin == 5) momtmp.SetZ(momtmp.Z() + delta);
105 Belle2::Helix helixPlusDelta(postmp, momtmp, charge, bfield);
107 jacobian_numerical(0, jin) = (helixPlusDelta.getD0() - helix.getD0()) / delta;
108 jacobian_numerical(1, jin) = (helixPlusDelta.getPhi0() - helix.getPhi0()) / delta;
109 jacobian_numerical(2, jin) = (helixPlusDelta.getOmega() - helix.getOmega()) / delta;
110 jacobian_numerical(3, jin) = (helixPlusDelta.getZ0() - helix.getZ0()) / delta;
111 jacobian_numerical(4, jin) = (helixPlusDelta.getTanLambda() - helix.getTanLambda()) / delta;
116 for (
int row = 0; row < 5; ++row) {
117 for (
int col = 0; col < 6; ++col) {
118 const double num = jacobian_numerical(row, col);
119 const double ana = jacobian_analytical(row, col);
120 const double res = std::abs(num - ana);
121 EXPECT_TRUE(res < eps) <<
"row " << row <<
" col " << col <<
" num - ana " << res <<
" num " << num <<
" ana " << ana;
static ROOT::Math::XYZVector getFieldInTesla(const ROOT::Math::XYZVector &pos)
return the magnetic field at a given position in Tesla.
static const double speedOfLight
[cm/ns]
static void getJacobianToCartesianFrameworkHelix(Eigen::Matrix< double, 5, 6 > &jacobian, const double x, const double y, const double z, const double px, const double py, const double pz, const double bfield, const double charge)
get the jacobian dh={helix pars}/dx={x,y,z,px,py,pz} for the implementation of the framework helix.
double charge(int pdgCode)
Returns electric charge of a particle with given pdg code.