Belle II Software light-2406-ragdoll
test_HelixJacobian.h
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8#pragma once
9
10#include <Eigen/Core>
11#include <gtest/gtest.h>
12
13#include <framework/gearbox/Const.h>
14
15#include <analysis/VertexFitting/TreeFitter/HelixUtils.h>
16#include <framework/geometry/BFieldManager.h>
17#include <framework/gearbox/Unit.h>
18
19namespace {
20
22 class TreeFitterHelixJacobianTest : public ::testing::Test {
23
24 };
25
27 TEST_F(TreeFitterHelixJacobianTest, Helix)
28 {
29 const double eps = 1E-6;
30 //random start values to ensure x and p vector are not orthogonal/parallel
31 const double px = 0.123214;
32 const double py = 1.543;
33 const double pz = -2.876;
34 const double x = 3.765346;
35 const double y = -2.56756;
36 const double z = 1.678;
37
38 const int charge = -1;
39 const double bfield = Belle2::BFieldManager::getFieldInTesla(ROOT::Math::XYZVector(0, 0, 0)).Z();
40 const double alpha = 1.0 / (bfield * Belle2::Const::speedOfLight) * 1E4;
41 const double aq = charge / alpha;
42
43 Belle2::Helix helix = Belle2::Helix(ROOT::Math::XYZVector(x, y, z), ROOT::Math::XYZVector(px, py, pz), charge, bfield);
44
45 const double pt = std::sqrt(px * px + py * py);
46 const double omega = aq / pt;
47
48 const double phi = atan(py / px);
49 const double cosPhi = std::cos(phi);
50 const double sinPhi = std::sin(phi);
51
52 const double para = -x * cosPhi - y * sinPhi;
53 const double ortho = -y * cosPhi + x * sinPhi;
54 const double R2 = para * para + ortho * ortho;
55 const double A = 2 * ortho + omega * R2;
56 const double U = std::sqrt(1 + omega * A);
57 const double d0 = A / (1 + U);
58 const double l = 1 / omega * atan((omega * para) / (1 + omega * ortho));
59 const double phi0 = phi + atan((omega * para) / (1 + omega * ortho));
60 const double tanLambda = pz / pt;
61 const double z0 = z + tanLambda * l;
62
63 std::vector<double> h = {d0, phi0, omega, z0, tanLambda};
64 std::vector<double> h_framework = {helix.getD0(), helix.getPhi0(), helix.getOmega(), helix.getZ0(), helix.getTanLambda()};
65 for (int row = 0; row < 5; ++row) {
66 double res = h[row] - h_framework[row];
67 EXPECT_TRUE(res < eps) << "row " << row << " num - ana " << res << " framework " << h_framework[row] << " mine " << h[row];
68 }
69 }
70
72 TEST_F(TreeFitterHelixJacobianTest, Parameters)
73 {
74 const double delta = 1e-6;
75 const double eps = 1e-5;
76
77 Eigen::Matrix<double, 5, 6> jacobian_numerical = Eigen::Matrix<double, 5, 6>::Zero(5, 6);
78 Eigen::Matrix<double, 5, 6> jacobian_analytical = Eigen::Matrix<double, 5, 6>::Zero(5, 6);
79
80 ROOT::Math::XYZVector postmp;
81 ROOT::Math::XYZVector momtmp;
82
83 //random start values to ensure x and p vector are not orthogonal/parallel
84 const double px = 0.523214;
85 const double py = -1.543;
86 const double pz = -2.876;
87 const double x = -3.765346;
88 const double y = -2.56756;
89 const double z = 5.678;
90
91 const Eigen::Matrix<double, 1, 6> positionAndMom_ = (Eigen::Matrix<double, 1, 6>() << x, y, z, px, py, pz).finished();
92 const int charge = -1;
93 const double bfield = Belle2::BFieldManager::getFieldInTesla(ROOT::Math::XYZVector(0, 0, 0)).Z();
94
95 Belle2::Helix helix = Belle2::Helix(ROOT::Math::XYZVector(x, y, z), ROOT::Math::XYZVector(px, py, pz), charge, bfield);
96
97 for (int jin = 0; jin < 6; ++jin) {
98 postmp.SetCoordinates(positionAndMom_(0), positionAndMom_(1), positionAndMom_(2));
99 momtmp.SetCoordinates(positionAndMom_(3), positionAndMom_(4), positionAndMom_(5));
100 if (jin == 0) postmp.SetX(postmp.X() + delta);
101 if (jin == 1) postmp.SetY(postmp.Y() + delta);
102 if (jin == 2) postmp.SetZ(postmp.Z() + delta);
103 if (jin == 3) momtmp.SetX(momtmp.X() + delta);
104 if (jin == 4) momtmp.SetY(momtmp.Y() + delta);
105 if (jin == 5) momtmp.SetZ(momtmp.Z() + delta);
106
107 Belle2::Helix helixPlusDelta(postmp, momtmp, charge, bfield);
108
109 jacobian_numerical(0, jin) = (helixPlusDelta.getD0() - helix.getD0()) / delta;
110 jacobian_numerical(1, jin) = (helixPlusDelta.getPhi0() - helix.getPhi0()) / delta;
111 jacobian_numerical(2, jin) = (helixPlusDelta.getOmega() - helix.getOmega()) / delta;
112 jacobian_numerical(3, jin) = (helixPlusDelta.getZ0() - helix.getZ0()) / delta;
113 jacobian_numerical(4, jin) = (helixPlusDelta.getTanLambda() - helix.getTanLambda()) / delta;
114 }
115
116 TreeFitter::HelixUtils::getJacobianToCartesianFrameworkHelix(jacobian_analytical, x, y, z, px, py, pz, bfield, charge);
117
118 for (int row = 0; row < 5; ++row) {
119 for (int col = 0; col < 6; ++col) {
120 const double num = jacobian_numerical(row, col);
121 const double ana = jacobian_analytical(row, col);
122 const double res = std::abs(num - ana);
123 EXPECT_TRUE(res < eps) << "row " << row << " col " << col << " num - ana " << res << " num " << num << " ana " << ana;
124 }
125 }
126 }
127
128}
static ROOT::Math::XYZVector getFieldInTesla(const ROOT::Math::XYZVector &pos)
return the magnetic field at a given position in Tesla.
Definition: BFieldManager.h:61
static const double speedOfLight
[cm/ns]
Definition: Const.h:695
This class represents an ideal helix in perigee parameterization.
Definition: Helix.h:59
double getOmega() const
Getter for omega, which is a signed curvature measure of the track.
Definition: Helix.h:387
double getD0() const
Getter for d0, which is the signed distance to the perigee in the r-phi plane.
Definition: Helix.h:372
double getTanLambda() const
Getter for tan lambda, which is the z over two dimensional arc length slope of the track.
Definition: Helix.h:393
double getZ0() const
Getter for z0, which is the z coordinate of the perigee.
Definition: Helix.h:390
double getPhi0() const
Getter for phi0, which is the azimuth angle of the transverse momentum at the perigee.
Definition: Helix.h:378
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.
Definition: HelixUtils.cc:407
double charge(int pdgCode)
Returns electric charge of a particle with given pdg code.
Definition: EvtPDLUtil.cc:44