Belle II Software development
HelixUtils.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * External Contributor: Wouter Hulsbergen *
5 * *
6 * See git log for contributors and copyright holders. *
7 * This file is licensed under LGPL-3.0, see LICENSE.md. *
8 **************************************************************************/
9
10#include <analysis/VertexFitting/TreeFitter/HelixUtils.h>
11
12#include <TMath.h>
13
14#include <framework/gearbox/Const.h>
15#include <framework/logging/Logger.h>
16#include <framework/dataobjects/Helix.h>
17
18#include <algorithm>
19#include <initializer_list>
20#include <vector>
21
22namespace TreeFitter {
23
24 void HelixUtils::vertexFromHelix(const Belle2::Helix& helix,
25 double L, double Bz,
26 ROOT::Math::XYZVector& position,
27 ROOT::Math::XYZVector& momentum, int& charge)
28 {
29 position = helix.getPositionAtArcLength2D(L);
30 momentum = helix.getMomentumAtArcLength2D(L, Bz);
31 charge = helix.getChargeSign();
32 }
33
34 void HelixUtils::helixFromVertex(const Eigen::Matrix<double, 1, 6>& positionAndMomentum,
35 int charge, double Bz,
36 Belle2::Helix& helix,
37 double& L,
38 Eigen::Matrix<double, 5, 6>& jacobian)
39 {
40
41 helix = Belle2::Helix(ROOT::Math::XYZVector(positionAndMomentum(0), positionAndMomentum(1), positionAndMomentum(2)),
42 ROOT::Math::XYZVector(positionAndMomentum(3), positionAndMomentum(4), positionAndMomentum(5)),
43 charge, Bz);
44
45 L = helix.getArcLength2DAtXY(positionAndMomentum(0),
46 positionAndMomentum(1));
47
48 const double alpha = helix.getAlpha(Bz);
49
50 //Copied from Belle2::UncertainHelix
51 // COMPLETELY WRONG SINCE IT ASSUMES IT'S IN THE.operator() PERIGEE,
52 // ONLY A PLACEHOLDER FOR NOW
53 // 1. Rotate to a system where phi0 = 0
54 Eigen::Matrix<double, 6, 6> jacobianRot = Eigen::Matrix<double, 6, 6>::Zero(6, 6);
55
56 const double px = positionAndMomentum(3);
57 const double py = positionAndMomentum(4);
58 const double pt = hypot(px, py);
59 const double cosPhi0 = px / pt;
60 const double sinPhi0 = py / pt;
61
62 // Passive rotation matrix by phi0:
63 jacobianRot(iX, iX) = cosPhi0;
64 jacobianRot(iX, iY) = sinPhi0;
65 jacobianRot(iY, iX) = -sinPhi0;
66 jacobianRot(iY, iY) = cosPhi0;
67 jacobianRot(iZ, iZ) = 1.0;
68
69 jacobianRot(iPx, iPx) = cosPhi0;
70 jacobianRot(iPx, iPy) = sinPhi0;
71 jacobianRot(iPy, iPx) = -sinPhi0;
72 jacobianRot(iPy, iPy) = cosPhi0;
73 jacobianRot(iPz, iPz) = 1.0;
74
75 // 2. Translate to perigee parameters on the position
76 const double pz = positionAndMomentum(5);
77 const double invPt = 1 / pt;
78 const double invPtSquared = invPt * invPt;
79 Eigen::Matrix<double, 5, 6> jacobianToHelixParameters = Eigen::Matrix<double, 5, 6>::Zero(5, 6);
80 jacobianToHelixParameters(iD0, iY) = -1;
81 jacobianToHelixParameters(iPhi0, iX) = charge * invPt / alpha;
82 jacobianToHelixParameters(iPhi0, iPy) = invPt;
83 jacobianToHelixParameters(iOmega, iPx) = -charge * invPtSquared / alpha;
84 jacobianToHelixParameters(iTanLambda, iPx) = - pz * invPtSquared;
85 jacobianToHelixParameters(iTanLambda, iPz) = invPt;
86 jacobianToHelixParameters(iZ0, iX) = - pz * invPt;
87 jacobianToHelixParameters(iZ0, iZ) = 1;
88 //
89 jacobian = jacobianToHelixParameters * jacobianRot;
90
91 }
92
93 std::string HelixUtils::helixParName(int i)
94 {
95 std::string rc ;
96 switch (i) {
97 case 1 : rc = "d0 : " ; break ;
98 case 2 : rc = "phi0 : " ; break ;
99 case 3 : rc = "omega : " ; break ;
100 case 4 : rc = "z0 : " ; break ;
101 case 5 : rc = "tandip: " ; break ;
102 case 6 : rc = "L : " ; break ;
103 }
104 return rc ;
105 }
106
107 std::string HelixUtils::vertexParName(int i)
108 {
109 std::string rc ;
110 switch (i) {
111 case 1 : rc = "x : " ; break ;
112 case 2 : rc = "y : " ; break ;
113 case 3 : rc = "z : " ; break ;
114 case 4 : rc = "px : " ; break ;
115 case 5 : rc = "py : " ; break ;
116 case 6 : rc = "pz : " ; break ;
117 }
118 return rc ;
119 }
120
121 void HelixUtils::printVertexPar(const ROOT::Math::XYZVector& position, const ROOT::Math::XYZVector& momentum, int charge)
122 {
123 B2INFO(vertexParName(1).c_str() << position.X());
124 B2INFO(vertexParName(2).c_str() << position.Y());
125 B2INFO(vertexParName(3).c_str() << position.Z());
126 B2INFO(vertexParName(4).c_str() << momentum.X());
127 B2INFO(vertexParName(5).c_str() << momentum.Y());
128 B2INFO(vertexParName(6).c_str() << momentum.Z());
129 B2INFO("charge: " << charge);
130
131 }
132
133 void HelixUtils::getHelixAndJacobianFromVertexNumerical(const Eigen::Matrix<double, 1, 6>& positionAndMom,
134 int charge, double Bz,
135 Belle2::Helix& helix,
136 // output parameter: Eigen's const operator() overload hides the
137 // assignments below from cppcheck
138 // cppcheck-suppress constParameterReference
139 Eigen::Matrix<double, 5, 6>& jacobian)
140 {
141
142 helix = Belle2::Helix(ROOT::Math::XYZVector(positionAndMom(0), positionAndMom(1), positionAndMom(2)),
143 ROOT::Math::XYZVector(positionAndMom(3), positionAndMom(4), positionAndMom(5)),
144 charge, Bz);
145
146 // numeric calculation of the jacobian
147 Belle2::Helix helixPlusDelta;
148
149 double delta = 1e-5;// this is arbitrary, only needs to be small
150
151 ROOT::Math::XYZVector postmp;
152 ROOT::Math::XYZVector momtmp;
153
154 for (int jin = 0; jin < 6; ++jin) {
155 postmp.SetCoordinates(positionAndMom(0), positionAndMom(1), positionAndMom(2));
156 momtmp.SetCoordinates(positionAndMom(3), positionAndMom(4), positionAndMom(5));
157 if (jin == 0) postmp.SetX(postmp.X() + delta);
158 if (jin == 1) postmp.SetY(postmp.Y() + delta);
159 if (jin == 2) postmp.SetZ(postmp.Z() + delta);
160 if (jin == 3) momtmp.SetX(momtmp.X() + delta);
161 if (jin == 4) momtmp.SetY(momtmp.Y() + delta);
162 if (jin == 5) momtmp.SetZ(momtmp.Z() + delta);
163
164 helixPlusDelta = Belle2::Helix(postmp, momtmp, charge, Bz);
165 jacobian(iD0, jin) = (helixPlusDelta.getD0() - helix.getD0()) / delta ;
166 jacobian(iPhi0, jin) = (helixPlusDelta.getPhi0() - helix.getPhi0()) / delta ;
167 jacobian(iOmega, jin) = (helixPlusDelta.getOmega() - helix.getOmega()) / delta ;
168 jacobian(iZ0, jin) = (helixPlusDelta.getZ0() - helix.getZ0()) / delta ;
169 jacobian(iTanLambda, jin) = (helixPlusDelta.getTanLambda() - helix.getTanLambda()) / delta ;
170
171 // jacobian[iArcLength2D][jin] = (LPlusDelta - L) / delta ;
172 }
173 }
174
175 void HelixUtils::getJacobianFromVertexNumerical(
176 const Eigen::Matrix<double, 1, 6>& positionAndMom,
177 int charge, double Bz,
178 const Belle2::Helix& helix,
179 // output parameter: Eigen's const operator() overload hides the assignments
180 // below from cppcheck
181 // cppcheck-suppress constParameterReference
182 Eigen::Matrix<double, 5, 6>& jacobian,
183 double delta
184 )
185 {
186 // numeric calculation of the jacobian
187 Belle2::Helix helixPlusDelta;
188
189 ROOT::Math::XYZVector postmp;
190 ROOT::Math::XYZVector momtmp;
191
192 for (int jin = 0; jin < 6; ++jin) {
193 postmp.SetCoordinates(positionAndMom(0), positionAndMom(1), positionAndMom(2));
194 momtmp.SetCoordinates(positionAndMom(3), positionAndMom(4), positionAndMom(5));
195 if (jin == 0) postmp.SetX(postmp.X() + delta);
196 if (jin == 1) postmp.SetY(postmp.Y() + delta);
197 if (jin == 2) postmp.SetZ(postmp.Z() + delta);
198 if (jin == 3) momtmp.SetX(momtmp.X() + delta);
199 if (jin == 4) momtmp.SetY(momtmp.Y() + delta);
200 if (jin == 5) momtmp.SetZ(momtmp.Z() + delta);
201
202 helixPlusDelta = Belle2::Helix(postmp, momtmp, charge, Bz);
203 jacobian(iD0, jin) = (helixPlusDelta.getD0() - helix.getD0()) / delta ;
204 jacobian(iPhi0, jin) = (helixPlusDelta.getPhi0() - helix.getPhi0()) / delta ;
205 jacobian(iOmega, jin) = (helixPlusDelta.getOmega() - helix.getOmega()) / delta ;
206 jacobian(iZ0, jin) = (helixPlusDelta.getZ0() - helix.getZ0()) / delta ;
207 jacobian(iTanLambda, jin) = (helixPlusDelta.getTanLambda() - helix.getTanLambda()) / delta ;
208 }
209
210 }
211
212 inline double sqr(double x) { return x * x ; }
213
214 double HelixUtils::phidomain(const double phi)
215 {
216 double rc = phi ;
217 if (phi < -TMath::Pi()) rc += TMath::TwoPi();
218 else if (phi > TMath::Pi()) rc -= TMath::TwoPi();
219 return rc ;
220 }
221
222 //POCA between two tracks
223 double HelixUtils::helixPoca(const Belle2::Helix& helix1,
224 const Belle2::Helix& helix2,
225 double& flt1, double& flt2,
226 Eigen::Vector3d& vertex, bool parallel)
227 {
228
229 const double d0_1 = helix1.getD0();
230 const double phi0_1 = helix1.getPhi0();
231 const double omega_1 = helix1.getOmega();
232
233 const double d0_2 = helix2.getD0();
234 const double phi0_2 = helix2.getPhi0();
235 const double omega_2 = helix2.getOmega();
236
237 // These radii have a sign, like omega (negative for negative charge)
238 const double r_1 = 1 / omega_1 ;
239 const double r_2 = 1 / omega_2 ;
240
241 // 1) First look at the transverse plane, where the helix projection is a circle
242 // Coordinates of the centers of the circles
243 const double x0_1 = (r_1 + d0_1) * sin(phi0_1) ;
244 const double y0_1 = -(r_1 + d0_1) * cos(phi0_1) ;
245
246 const double x0_2 = (r_2 + d0_2) * sin(phi0_2) ;
247 const double y0_2 = -(r_2 + d0_2) * cos(phi0_2) ;
248
249 // Vector that goes from center1 to center2
250 const double deltax = x0_2 - x0_1 ;
251 const double deltay = y0_2 - y0_1 ;
252
253 // Intersections of the circles, can be at most two
254 double phi1[2] ;
255 double phi2[2] ;
256 int nsolutions = 1;
257
258 // The phi of the delta vector.
259 const double phi = - atan2(deltax, deltay) ;
260 const double phinot = phi > 0 ? phi - TMath::Pi() : phi + TMath::Pi() ;
261 phi1[0] = r_1 < 0 ? phi : phinot ;
262 phi2[0] = r_2 > 0 ? phi : phinot ;
263
264 // These radii do NOT have a sign instead
265 const double R1 = fabs(r_1) ;
266 const double R2 = fabs(r_2) ;
267 const double Rmin = R1 < R2 ? R1 : R2 ;
268 const double Rmax = R1 > R2 ? R1 : R2 ;
269 const double dX = hypot(deltax, deltay) ;
270
271 if (!parallel && dX + Rmin > Rmax && dX < R1 + R2) {
272 // Circles intersect in two points
273 nsolutions = 2 ;
274
275 // This is just the law of cosines
276 const double ddphi1 = acos((dX * dX - R2 * R2 + R1 * R1) / (2.*dX * R1)) ;
277 phi1[1] = phidomain(phi1[0] + ddphi1) ;
278 phi1[0] = phidomain(phi1[0] - ddphi1) ;
279
280 const double ddphi2 = acos((dX * dX - R1 * R1 + R2 * R2) / (2.*dX * R2)) ;
281 phi2[1] = phidomain(phi2[0] - ddphi2) ;
282 phi2[0] = phidomain(phi2[0] + ddphi2) ;
283
284 } else if (dX < Rmax) {
285 // Tangent or non-intersecting circles, one inside the other (only one POCA)
286 // the two branches are symmetric; keep both assignments explicit
287 if (R1 > R2) phi2[0] = r_2 < 0 ? phi : phinot ;
288 // cppcheck-suppress redundantAssignment
289 else phi1[0] = r_1 < 0 ? phi : phinot ;
290 }
291 // else: tangent or non-intersecting circles, outside of each other (only one POCA)
292 // what we saved in phi1 and phi2 gives already the correct solution
293
294 // Intersections of the circles (cartesian)
295 double x1[2], y1[2], x2[2], y2[2];
296 for (int i = 0; i < nsolutions; i++) {
297 x1[i] = r_1 * sin(phi1[i]) + x0_1 ;
298 y1[i] = -r_1 * cos(phi1[i]) + y0_1 ;
299 x2[i] = r_2 * sin(phi2[i]) + x0_2 ;
300 y2[i] = -r_2 * cos(phi2[i]) + y0_2 ;
301 }
302
303 // 2) Find the best solution for z by running multiples of 2pi from the xy intersection(s)
304 double z1(0), z2(0);
305 bool first = true;
306 int ibest = 0;
307 const int nturnsmax = 10; // Max number of turns we try backwards and forwards
308
309 // Loop on all xy-plane solutions
310 for (int i = 0; i < nsolutions; ++i) {
311 const double l1 = helix1.getArcLength2DAtXY(x1[i], y1[i]);
312 const double l2 = helix2.getArcLength2DAtXY(x2[i], y2[i]);
313
314 // Loop on helix1 turns, save corresponding z positions
315 std::vector<double> z1s;
316 for (int n1 = 0; n1 <= nturnsmax; ++n1) {
317 bool added = false;
318 // Try forwards and backwards
319 for (int sn1 : {n1, -n1}) {
320 const double tmpz1 = helix1.getPositionAtArcLength2D(l1 + sn1 * TMath::TwoPi() / omega_1).Z();
321 if (sn1 == 0 || (-82 <= tmpz1 && tmpz1 <= 158)) {
322 // Only keep the 0th turn and those inside CDC volume
323 z1s.push_back(tmpz1);
324 added = true;
325 }
326 if (sn1 == 0)
327 break; // Do not store 0th turn twice
328 }
329 // If we did not add any point we are already outside CDC volume both backwards and forwards
330 if (!added)
331 break;
332 }
333
334 // Loop on helix2 turns, find closest approach to one of helix1 points
335 for (int n2 = 0; n2 <= nturnsmax; ++n2) {
336 bool tried = false;
337 // Try forwards and backwards
338 for (int sn2 : {n2, -n2}) {
339 const double tmpz2 = helix2.getPositionAtArcLength2D(l2 + sn2 * TMath::TwoPi() / omega_2).Z();
340 if (sn2 == 0 || (-82 <= tmpz2 && tmpz2 <= 158)) {
341 // Only keep the 0th turn and those inside CDC volume
342 tried = true;
343 // Find the tmpz1 closest to tmpz2
344 const auto i1best = std::min_element(
345 z1s.cbegin(), z1s.cend(), [&tmpz2](const double & z1a, const double & z1b) {
346 return fabs(z1a - tmpz2) < fabs(z1b - tmpz2);
347 });
348 const double tmpz1 = *i1best;
349 // Keep the solution where the z distance of closest approach is minimum
350 if (first || fabs(tmpz1 - tmpz2) < fabs(z1 - z2)) {
351 ibest = i;
352 first = false;
353 z1 = tmpz1;
354 z2 = tmpz2;
355 flt1 = l1;
356 flt2 = l2;
357 }
358 }
359 if (n2 == 0)
360 break; // Do not try 0th turn twice
361 }
362 // If we did not try any point we are already outside CDC volume both backwards and forwards
363 if (!tried)
364 break;
365 }
366 }
367
368 vertex.x() = 0.5 * (x1[ibest] + x2[ibest]);
369 vertex.y() = 0.5 * (y1[ibest] + y2[ibest]);
370 vertex.z() = 0.5 * (z1 + z2);
371
372 return std::hypot(x2[ibest] - x1[ibest], y2[ibest] - y1[ibest], z2 - z1);
373 }
374
375 //POCA between a track and a point
376 double HelixUtils::helixPoca(const Belle2::Helix& helix,
377 const ROOT::Math::XYZVector& point,
378 double& flt)
379 {
380 const double d0 = helix.getD0();
381 const double phi0 = helix.getPhi0();
382 const double omega = helix.getOmega();
383 const double z0 = helix.getZ0();
384 const double tandip = helix.getTanLambda();
385 const double cosdip = cos(atan(tandip)) ; // can do that faster
386
387 const double r = 1 / omega ;
388
389 const double x0 = - (r + d0) * sin(phi0) ;
390 const double y0 = (r + d0) * cos(phi0) ;
391
392 const double deltax = x0 - point.X() ;
393 const double deltay = y0 - point.Y() ;
394
395 const double pi = TMath::Pi();
396 double phi = - atan2(deltax, deltay) ;
397 if (r < 0) phi = phi > 0 ? phi - pi : phi + pi ;
398
399 // find the best solution for z by running multiples of 2_pi
400 const double x = r * sin(phi) + x0 ;
401 const double y = -r * cos(phi) + y0 ;
402 double z(0) ;
403 bool first(true) ;
404 const int ncirc(2) ;
405 const double dphi = phidomain(phi - phi0) ;
406 for (int n = 1 - ncirc; n <= 1 + ncirc ; ++n) {
407 const double l = (dphi + n * TMath::TwoPi()) / omega ;
408 const double tmpz = (z0 + l * tandip) ;
409 if (first || fabs(tmpz - point.Z()) < fabs(z - point.Z())) {
410 first = false ;
411 z = tmpz ;
412 flt = l / cosdip ;
413 }
414 }
415 return sqrt(sqr(x - point.X()) + sqr(y - point.Y()) + sqr(z - point.Z())) ;
416 }
417
418 // cppcheck-suppress constParameterReference ; output parameter, see above
419 void HelixUtils::getJacobianToCartesianFrameworkHelix(Eigen::Matrix<double, 5, 6>& jacobian,
420 const double x,
421 const double y,
422 const double z __attribute__((unused)),
423 const double px,
424 const double py,
425 const double pz,
426 const double bfield,
427 const double charge
428 )
429
430 {
431 const double alpha = 1.0 / (bfield * Belle2::Const::speedOfLight) * 1E4;
432 const double aq = charge / alpha;
433
434 const double pt = std::hypot(px, py);
435 const double pt2 = pt * pt;
436 const double pt3 = pt2 * pt;
437 const double aq2 = aq * aq;
438
439 const double x2 = x * x;
440 const double y2 = y * y;
441 const double r = x2 + y2;
442
443 const double px2 = px * px;
444 const double py2 = py * py;
445
446 const double px0 = px - aq * y;
447 const double py0 = py + aq * x;
448
449 const double pt02 = px0 * px0 + py0 * py0;
450 const double pt0 = std::sqrt(pt02);
451 double sqrt13 = pt0 / pt;
452
453 // D d0 / Dx_i
454 jacobian(0, 0) = py0 / pt0;
455 jacobian(0, 1) = -px0 / pt0;
456 jacobian(0, 2) = 0;
457 jacobian(0, 3) = (-(y * (aq2 * r + 2 * aq * py * x + 2 * py2 * (1 + sqrt13))) - px * (2 * py * x * (1 + sqrt13) + aq * (y2 *
458 (-1 + sqrt13) + x2 * (1 + sqrt13)))) /
459 (pt2 * pt0 * (1 + sqrt13) * (1 + sqrt13));
460
461 jacobian(0, 4) = (2 * px2 * x * (1 + sqrt13) + 2 * px * y * (py - aq * x + py * sqrt13) + aq * (aq * r * x - py * (x2 *
462 (-1 + sqrt13) + y2 * (1 + sqrt13)))) /
463 (pt2 * pt0 * (1 + sqrt13) * (1 + sqrt13));
464 jacobian(0, 5) = 0;
465
466 // D phi0 / Dx_i0;
467 jacobian(1, 0) = aq * px0 / pt02;
468 jacobian(1, 1) = aq * py0 / pt02;
469 jacobian(1, 2) = 0;
470 jacobian(1, 3) = -py0 / pt02;
471 jacobian(1, 4) = px0 / pt02;
472 jacobian(1, 5) = 0;
473
474 // D omega / Dx_i
475 jacobian(2, 0) = 0;
476 jacobian(2, 1) = 0;
477 jacobian(2, 2) = 0;
478 jacobian(2, 3) = - aq * px / pt3;
479 jacobian(2, 4) = - aq * py / pt3;
480 jacobian(2, 5) = 0;
481
482 // D z0 / Dx_i
483 jacobian(3, 0) = -pz * px0 / pt02;
484 jacobian(3, 1) = -pz * py0 / pt02;
485 jacobian(3, 2) = 1;
486 jacobian(3, 3) = (pz * (px2 * x - py * (aq * r + py * x) + 2 * px * py * y)) / (pt2 * pt02);
487 jacobian(3, 4) = (pz * (px * (aq * r + 2 * py * x) - px2 * y + py2 * y)) / (pt2 * pt02);
488 jacobian(3, 5) = std::atan2(-(aq * (px * x + py * y)), (px2 + py * py0 - aq * px * y)) / aq; //pt on num. and denom cancels.
489
490 // D tan lambda / Dx_i
491 jacobian(4, 0) = 0;
492 jacobian(4, 1) = 0;
493 jacobian(4, 2) = 0;
494 jacobian(4, 3) = -pz * px / pt3;
495 jacobian(4, 4) = -pz * py / pt3;
496 jacobian(4, 5) = 1. / pt;
497 }
498
499}
Belle2::CDC::Helix
Helix parameter class.
Definition Helix.h:48
Belle2::Const::speedOfLight
static const double speedOfLight
[cm/ns]
Definition Const.h:696
TreeFitter::HelixUtils::vertexParName
static std::string vertexParName(int i)
map of the vertex parameters by list index
Definition HelixUtils.cc:107
TreeFitter::HelixUtils::printVertexPar
static void printVertexPar(const ROOT::Math::XYZVector &position, const ROOT::Math::XYZVector &momentum, int charge)
Print the vertex parameters.
Definition HelixUtils.cc:121
TreeFitter::HelixUtils::helixFromVertex
static void helixFromVertex(const Eigen::Matrix< double, 1, 6 > &positionAndMomentum, int charge, double Bz, Belle2::Helix &helix, double &L, Eigen::Matrix< double, 5, 6 > &jacobian)
vertex --> helix
Definition HelixUtils.cc:34
TreeFitter::HelixUtils::getHelixAndJacobianFromVertexNumerical
static void getHelixAndJacobianFromVertexNumerical(const Eigen::Matrix< double, 1, 6 > &positionAndMom, int charge, double Bz, Belle2::Helix &helix, Eigen::Matrix< double, 5, 6 > &jacobian)
get helix and jacobian from a vertex
Definition HelixUtils.cc:133
TreeFitter::HelixUtils::vertexFromHelix
static void vertexFromHelix(const Belle2::Helix &helix, double L, double Bz, ROOT::Math::XYZVector &position, ROOT::Math::XYZVector &momentum, int &charge)
helix --> vertex
Definition HelixUtils.cc:24
TreeFitter::HelixUtils::getJacobianFromVertexNumerical
static void getJacobianFromVertexNumerical(const Eigen::Matrix< double, 1, 6 > &positionAndMom, int charge, double Bz, const Belle2::Helix &helix, Eigen::Matrix< double, 5, 6 > &jacobian, double delta=1e-5)
get jacobian from a vertex
Definition HelixUtils.cc:175
TreeFitter::HelixUtils::helixParName
static std::string helixParName(int i)
map of the helix parameters by list index
Definition HelixUtils.cc:93
TreeFitter::HelixUtils::phidomain
static double phidomain(const double phi)
the domain of phi
Definition HelixUtils.cc:214
TreeFitter::HelixUtils::getJacobianToCartesianFrameworkHelix
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:419
TreeFitter::HelixUtils::helixPoca
static double helixPoca(const Belle2::Helix &helix1, const Belle2::Helix &helix2, double &flt1, double &flt2, Eigen::Vector3d &vertex, bool parallel=false)
POCA between two tracks.
Definition HelixUtils.cc:223