Belle II Software development
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 <gtest/gtest.h>
11
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>
16
17namespace {
18
20 class TreeFitterHelixJacobianTest : public ::testing::Test {
21
22 };
23
25 TEST_F(TreeFitterHelixJacobianTest, Helix)
26 {
27 const double eps = 1E-6;
28 //random start values to ensure x and p vector are not orthogonal/parallel
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;
35
36 const int charge = -1;
37 const double bfield = Belle2::BFieldManager::getFieldInTesla(ROOT::Math::XYZVector(0, 0, 0)).Z();
38 const double alpha = 1.0 / (bfield * Belle2::Const::speedOfLight) * 1E4;
39 const double aq = charge / alpha;
40
41 Belle2::Helix helix = Belle2::Helix(ROOT::Math::XYZVector(x, y, z), ROOT::Math::XYZVector(px, py, pz), charge, bfield);
42
43 const double pt = std::sqrt(px * px + py * py);
44 const double omega = aq / pt;
45
46 const double phi = atan(py / px);
47 const double cosPhi = std::cos(phi);
48 const double sinPhi = std::sin(phi);
49
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;
60
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];
66 }
67 }
68
70 TEST_F(TreeFitterHelixJacobianTest, Parameters)
71 {
72 const double delta = 1e-6;
73 const double eps = 1e-5;
74
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);
77
78 ROOT::Math::XYZVector postmp;
79 ROOT::Math::XYZVector momtmp;
80
81 //random start values to ensure x and p vector are not orthogonal/parallel
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;
88
89 const Eigen::Matrix<double, 1, 6> positionAndMom_ = (Eigen::Matrix<double, 1, 6>() << x, y, z, px, py, pz).finished();
90 const int charge = -1;
91 const double bfield = Belle2::BFieldManager::getFieldInTesla(ROOT::Math::XYZVector(0, 0, 0)).Z();
92
93 Belle2::Helix helix = Belle2::Helix(ROOT::Math::XYZVector(x, y, z), ROOT::Math::XYZVector(px, py, pz), charge, bfield);
94
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);
104
105 Belle2::Helix helixPlusDelta(postmp, momtmp, charge, bfield);
106
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;
112 }
113
114 TreeFitter::HelixUtils::getJacobianToCartesianFrameworkHelix(jacobian_analytical, x, y, z, px, py, pz, bfield, charge);
115
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;
122 }
123 }
124 }
125
126}
R E
internal precision of FFTW codelets
static ROOT::Math::XYZVector getFieldInTesla(const ROOT::Math::XYZVector &pos)
return the magnetic field at a given position in Tesla.
Helix parameter class.
Definition Helix.h:48
static const double speedOfLight
[cm/ns]
Definition Const.h:695
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.
Definition EvtPDLUtil.cc:44