Belle II Software light-2406-ragdoll
HelixUtils Class Reference

utility for helix<->x,p conversions More...

#include <HelixUtils.h>

Public Types

enum  VertexCoor {
  iX = 0 ,
  iY ,
  iZ ,
  iPx ,
  iPy ,
  iPz
}
 Parameters of the vertex
More...
 
enum  HelixCoor {
  iD0 = 0 ,
  iPhi0 ,
  iOmega ,
  iZ0 ,
  iTanLambda ,
  iArcLength2D
}
 Parameters of the helix. More...
 

Static Public Member Functions

static void helixFromVertexNumerical (const Eigen::Matrix< double, 3, 1 > &position, const Eigen::Matrix< double, 3, 1 > &momentum, int charge, double Bz, Belle2::Helix &helix, double &flt, Eigen::Matrix< double, 5, 6 > &jacobian)
 get helix from a vertex
 
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.
 
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
 
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
 
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
 
static void vertexFromHelix (const Belle2::Helix &helix, double L, double Bz, ROOT::Math::XYZVector &position, ROOT::Math::XYZVector &momentum, int &charge)
 helix --> vertex
 
static std::string helixParName (int i)
 map of the helix parameters by list index
 
static std::string vertexParName (int i)
 map of the vertex parameters by list index
 
static void printVertexPar (const Belle2::B2Vector3D &position, const Belle2::B2Vector3D &momentum, int charge)
 Print the vertex parameters.
 
static double helixPoca (const Belle2::Helix &helix1, const Belle2::Helix &helix2, double &flt1, double &flt2, Belle2::B2Vector3D &vertex, bool parallel=false)
 POCA between two tracks.
 
static double helixPoca (const Belle2::Helix &helix, const Belle2::B2Vector3D &point, double &flt)
 POCA between a track and a point.
 
static double phidomain (const double phi)
 the domain of phi
 

Detailed Description

utility for helix<->x,p conversions

Definition at line 20 of file HelixUtils.h.

Member Enumeration Documentation

◆ HelixCoor

enum HelixCoor

Parameters of the helix.

Definition at line 28 of file HelixUtils.h.

28{iD0 = 0, iPhi0, iOmega, iZ0, iTanLambda, iArcLength2D} ;

◆ VertexCoor

enum VertexCoor

Parameters of the vertex

Definition at line 25 of file HelixUtils.h.

25{iX = 0, iY, iZ, iPx, iPy, iPz} ;

Member Function Documentation

◆ getHelixAndJacobianFromVertexNumerical()

void getHelixAndJacobianFromVertexNumerical ( const Eigen::Matrix< double, 1, 6 > &  positionAndMom,
int  charge,
double  Bz,
Belle2::Helix helix,
Eigen::Matrix< double, 5, 6 > &  jacobian 
)
static

get helix and jacobian from a vertex

Definition at line 130 of file HelixUtils.cc.

134 {
135
136 helix = Belle2::Helix(ROOT::Math::XYZVector(positionAndMom(0), positionAndMom(1), positionAndMom(2)),
137 ROOT::Math::XYZVector(positionAndMom(3), positionAndMom(4), positionAndMom(5)),
138 charge, Bz);
139
140 // numeric calculation of the jacobian
141 Belle2::Helix helixPlusDelta;
142
143 double delta = 1e-5;// this is arbitrary, only needs to be small
144
145 ROOT::Math::XYZVector postmp;
146 ROOT::Math::XYZVector momtmp;
147
148 for (int jin = 0; jin < 6; ++jin) {
149 postmp.SetCoordinates(positionAndMom(0), positionAndMom(1), positionAndMom(2));
150 momtmp.SetCoordinates(positionAndMom(3), positionAndMom(4), positionAndMom(5));
151 if (jin == 0) postmp.SetX(postmp.X() + delta);
152 if (jin == 1) postmp.SetY(postmp.Y() + delta);
153 if (jin == 2) postmp.SetZ(postmp.Z() + delta);
154 if (jin == 3) momtmp.SetX(momtmp.X() + delta);
155 if (jin == 4) momtmp.SetY(momtmp.Y() + delta);
156 if (jin == 5) momtmp.SetZ(momtmp.Z() + delta);
157
158 helixPlusDelta = Belle2::Helix(postmp, momtmp, charge, Bz);
159 jacobian(iD0, jin) = (helixPlusDelta.getD0() - helix.getD0()) / delta ;
160 jacobian(iPhi0, jin) = (helixPlusDelta.getPhi0() - helix.getPhi0()) / delta ;
161 jacobian(iOmega, jin) = (helixPlusDelta.getOmega() - helix.getOmega()) / delta ;
162 jacobian(iZ0, jin) = (helixPlusDelta.getZ0() - helix.getZ0()) / delta ;
163 jacobian(iTanLambda, jin) = (helixPlusDelta.getTanLambda() - helix.getTanLambda()) / delta ;
164
165 // jacobian[iArcLength2D][jin] = (LPlusDelta - L) / delta ;
166 }
167 }
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

◆ getJacobianFromVertexNumerical()

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 
)
static

get jacobian from a vertex

Definition at line 169 of file HelixUtils.cc.

176 {
177 // numeric calculation of the jacobian
178 Belle2::Helix helixPlusDelta;
179
180 ROOT::Math::XYZVector postmp;
181 ROOT::Math::XYZVector momtmp;
182
183 for (int jin = 0; jin < 6; ++jin) {
184 postmp.SetCoordinates(positionAndMom(0), positionAndMom(1), positionAndMom(2));
185 momtmp.SetCoordinates(positionAndMom(3), positionAndMom(4), positionAndMom(5));
186 if (jin == 0) postmp.SetX(postmp.X() + delta);
187 if (jin == 1) postmp.SetY(postmp.Y() + delta);
188 if (jin == 2) postmp.SetZ(postmp.Z() + delta);
189 if (jin == 3) momtmp.SetX(momtmp.X() + delta);
190 if (jin == 4) momtmp.SetY(momtmp.Y() + delta);
191 if (jin == 5) momtmp.SetZ(momtmp.Z() + delta);
192
193 helixPlusDelta = Belle2::Helix(postmp, momtmp, charge, Bz);
194 jacobian(iD0, jin) = (helixPlusDelta.getD0() - helix.getD0()) / delta ;
195 jacobian(iPhi0, jin) = (helixPlusDelta.getPhi0() - helix.getPhi0()) / delta ;
196 jacobian(iOmega, jin) = (helixPlusDelta.getOmega() - helix.getOmega()) / delta ;
197 jacobian(iZ0, jin) = (helixPlusDelta.getZ0() - helix.getZ0()) / delta ;
198 jacobian(iTanLambda, jin) = (helixPlusDelta.getTanLambda() - helix.getTanLambda()) / delta ;
199 }
200
201 }

◆ getJacobianToCartesianFrameworkHelix()

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 
)
static

get the jacobian dh={helix pars}/dx={x,y,z,px,py,pz} for the implementation of the framework helix.

WARNING only valid right after initialisation!

Definition at line 407 of file HelixUtils.cc.

418 {
419 const double alpha = 1.0 / (bfield * Belle2::Const::speedOfLight) * 1E4;
420 const double aq = charge / alpha;
421
422 const double pt = std::hypot(px, py);
423 const double pt2 = pt * pt;
424 const double pt3 = pt2 * pt;
425 const double aq2 = aq * aq;
426
427 const double x2 = x * x;
428 const double y2 = y * y;
429 const double r = x2 + y2;
430
431 const double px2 = px * px;
432 const double py2 = py * py;
433
434 const double px0 = px - aq * y;
435 const double py0 = py + aq * x;
436
437 const double pt02 = px0 * px0 + py0 * py0;
438 const double pt0 = std::sqrt(pt02);
439 double sqrt13 = pt0 / pt;
440
441 // D d0 / Dx_i
442 jacobian(0, 0) = py0 / pt0;
443 jacobian(0, 1) = -px0 / pt0;
444 jacobian(0, 2) = 0;
445 jacobian(0, 3) = (-(y * (aq2 * r + 2 * aq * py * x + 2 * py2 * (1 + sqrt13))) - px * (2 * py * x * (1 + sqrt13) + aq * (y2 *
446 (-1 + sqrt13) + x2 * (1 + sqrt13)))) /
447 (pt2 * pt0 * (1 + sqrt13) * (1 + sqrt13));
448
449 jacobian(0, 4) = (2 * px2 * x * (1 + sqrt13) + 2 * px * y * (py - aq * x + py * sqrt13) + aq * (aq * r * x - py * (x2 *
450 (-1 + sqrt13) + y2 * (1 + sqrt13)))) /
451 (pt2 * pt0 * (1 + sqrt13) * (1 + sqrt13));
452 jacobian(0, 5) = 0;
453
454 // D phi0 / Dx_i0;
455 jacobian(1, 0) = aq * px0 / pt02;
456 jacobian(1, 1) = aq * py0 / pt02;
457 jacobian(1, 2) = 0;
458 jacobian(1, 3) = -py0 / pt02;
459 jacobian(1, 4) = px0 / pt02;
460 jacobian(1, 5) = 0;
461
462 // D omega / Dx_i
463 jacobian(2, 0) = 0;
464 jacobian(2, 1) = 0;
465 jacobian(2, 2) = 0;
466 jacobian(2, 3) = - aq * px / pt3;
467 jacobian(2, 4) = - aq * py / pt3;
468 jacobian(2, 5) = 0;
469
470 // D z0 / Dx_i
471 jacobian(3, 0) = -pz * px0 / pt02;
472 jacobian(3, 1) = -pz * py0 / pt02;
473 jacobian(3, 2) = 1;
474 jacobian(3, 3) = (pz * (px2 * x - py * (aq * r + py * x) + 2 * px * py * y)) / (pt2 * pt02);
475 jacobian(3, 4) = (pz * (px * (aq * r + 2 * py * x) - px2 * y + py2 * y)) / (pt2 * pt02);
476 jacobian(3, 5) = std::atan2(-(aq * (px * x + py * y)), (px2 + py * py0 - aq * px * y)) / aq; //pt on num. and denom cancels.
477
478 // D tan lambda / Dx_i
479 jacobian(4, 0) = 0;
480 jacobian(4, 1) = 0;
481 jacobian(4, 2) = 0;
482 jacobian(4, 3) = -pz * px / pt3;
483 jacobian(4, 4) = -pz * py / pt3;
484 jacobian(4, 5) = 1. / pt;
485 }
static const double speedOfLight
[cm/ns]
Definition: Const.h:695
double charge(int pdgCode)
Returns electric charge of a particle with given pdg code.
Definition: EvtPDLUtil.cc:44

◆ helixFromVertex()

void helixFromVertex ( const Eigen::Matrix< double, 1, 6 > &  positionAndMomentum,
int  charge,
double  Bz,
Belle2::Helix helix,
double &  L,
Eigen::Matrix< double, 5, 6 > &  jacobian 
)
static

vertex --> helix

Definition at line 33 of file HelixUtils.cc.

38 {
39
40 helix = Belle2::Helix(ROOT::Math::XYZVector(positionAndMomentum(0), positionAndMomentum(1), positionAndMomentum(2)),
41 ROOT::Math::XYZVector(positionAndMomentum(3), positionAndMomentum(4), positionAndMomentum(5)),
42 charge, Bz);
43
44 L = helix.getArcLength2DAtXY(positionAndMomentum(0),
45 positionAndMomentum(1));
46
47 const double alpha = helix.getAlpha(Bz);
48
49 //Copied from Belle2::UncertainHelix
50 // COMPLETELY WRONG SINCE IT ASSUMES IT'S IN THE.operator() PERIGEE,
51 // ONLY A PLACEHOLDER FOR NOW
52 // 1. Rotate to a system where phi0 = 0
53 Eigen::Matrix<double, 6, 6> jacobianRot = Eigen::Matrix<double, 6, 6>::Zero(6, 6);
54
55 const double px = positionAndMomentum(3);
56 const double py = positionAndMomentum(4);
57 const double pt = hypot(px, py);
58 const double cosPhi0 = px / pt;
59 const double sinPhi0 = py / pt;
60
61 // Passive rotation matrix by phi0:
62 jacobianRot(iX, iX) = cosPhi0;
63 jacobianRot(iX, iY) = sinPhi0;
64 jacobianRot(iY, iX) = -sinPhi0;
65 jacobianRot(iY, iY) = cosPhi0;
66 jacobianRot(iZ, iZ) = 1.0;
67
68 jacobianRot(iPx, iPx) = cosPhi0;
69 jacobianRot(iPx, iPy) = sinPhi0;
70 jacobianRot(iPy, iPx) = -sinPhi0;
71 jacobianRot(iPy, iPy) = cosPhi0;
72 jacobianRot(iPz, iPz) = 1.0;
73
74 // 2. Translate to perigee parameters on the position
75 const double pz = positionAndMomentum(5);
76 const double invPt = 1 / pt;
77 const double invPtSquared = invPt * invPt;
78 Eigen::Matrix<double, 5, 6> jacobianToHelixParameters = Eigen::Matrix<double, 5, 6>::Zero(5, 6);
79 jacobianToHelixParameters(iD0, iY) = -1;
80 jacobianToHelixParameters(iPhi0, iX) = charge * invPt / alpha;
81 jacobianToHelixParameters(iPhi0, iPy) = invPt;
82 jacobianToHelixParameters(iOmega, iPx) = -charge * invPtSquared / alpha;
83 jacobianToHelixParameters(iTanLambda, iPx) = - pz * invPtSquared;
84 jacobianToHelixParameters(iTanLambda, iPz) = invPt;
85 jacobianToHelixParameters(iZ0, iX) = - pz * invPt;
86 jacobianToHelixParameters(iZ0, iZ) = 1;
87 //
88 jacobian = jacobianToHelixParameters * jacobianRot;
89
90 }
double getArcLength2DAtXY(const double &x, const double &y) const
Calculates the two dimensional arc length at which the circle in the xy projection is closest to the ...
Definition: Helix.cc:141
static double getAlpha(const double bZ)
Calculates the alpha value for a given magnetic field in Tesla.
Definition: Helix.cc:109

◆ helixParName()

std::string helixParName ( int  i)
static

map of the helix parameters by list index

Definition at line 92 of file HelixUtils.cc.

93 {
94 std::string rc ;
95 switch (i) {
96 case 1 : rc = "d0 : " ; break ;
97 case 2 : rc = "phi0 : " ; break ;
98 case 3 : rc = "omega : " ; break ;
99 case 4 : rc = "z0 : " ; break ;
100 case 5 : rc = "tandip: " ; break ;
101 case 6 : rc = "L : " ; break ;
102 }
103 return rc ;
104 }

◆ helixPoca() [1/2]

double helixPoca ( const Belle2::Helix helix,
const Belle2::B2Vector3D point,
double &  flt 
)
static

POCA between a track and a point.

Definition at line 365 of file HelixUtils.cc.

368 {
369 const double d0 = helix.getD0();
370 const double phi0 = helix.getPhi0();
371 const double omega = helix.getOmega();
372 const double z0 = helix.getZ0();
373 const double tandip = helix.getTanLambda();
374 const double cosdip = cos(atan(tandip)) ; // can do that faster
375
376 const double r = 1 / omega ;
377
378 const double x0 = - (r + d0) * sin(phi0) ;
379 const double y0 = (r + d0) * cos(phi0) ;
380
381 const double deltax = x0 - point.X() ;
382 const double deltay = y0 - point.Y() ;
383
384 const double pi = TMath::Pi();
385 double phi = - atan2(deltax, deltay) ;
386 if (r < 0) phi = phi > 0 ? phi - pi : phi + pi ;
387
388 // find the best solution for z by running multiples of 2_pi
389 const double x = r * sin(phi) + x0 ;
390 const double y = -r * cos(phi) + y0 ;
391 double z(0) ;
392 bool first(true) ;
393 const int ncirc(2) ;
394 const double dphi = phidomain(phi - phi0) ;
395 for (int n = 1 - ncirc; n <= 1 + ncirc ; ++n) {
396 const double l = (dphi + n * TMath::TwoPi()) / omega ;
397 const double tmpz = (z0 + l * tandip) ;
398 if (first || fabs(tmpz - point.Z()) < fabs(z - point.Z())) {
399 first = false ;
400 z = tmpz ;
401 flt = l / cosdip ;
402 }
403 }
404 return sqrt(sqr(x - point.X()) + sqr(y - point.Y()) + sqr(z - point.Z())) ;
405 }
DataType Z() const
access variable Z (= .at(2) without boundary check)
Definition: B2Vector3.h:435
DataType X() const
access variable X (= .at(0) without boundary check)
Definition: B2Vector3.h:431
DataType Y() const
access variable Y (= .at(1) without boundary check)
Definition: B2Vector3.h:433
static double phidomain(const double phi)
the domain of phi
Definition: HelixUtils.cc:205

◆ helixPoca() [2/2]

double helixPoca ( const Belle2::Helix helix1,
const Belle2::Helix helix2,
double &  flt1,
double &  flt2,
Belle2::B2Vector3D vertex,
bool  parallel = false 
)
static

POCA between two tracks.

Definition at line 214 of file HelixUtils.cc.

218 {
219
220 const double d0_1 = helix1.getD0();
221 const double phi0_1 = helix1.getPhi0();
222 const double omega_1 = helix1.getOmega();
223
224 const double d0_2 = helix2.getD0();
225 const double phi0_2 = helix2.getPhi0();
226 const double omega_2 = helix2.getOmega();
227
228 // These radii have a sign, like omega (negative for negative charge)
229 const double r_1 = 1 / omega_1 ;
230 const double r_2 = 1 / omega_2 ;
231
232 // 1) First look at the transverse plane, where the helix projection is a circle
233 // Coordinates of the centers of the circles
234 const double x0_1 = (r_1 + d0_1) * sin(phi0_1) ;
235 const double y0_1 = -(r_1 + d0_1) * cos(phi0_1) ;
236
237 const double x0_2 = (r_2 + d0_2) * sin(phi0_2) ;
238 const double y0_2 = -(r_2 + d0_2) * cos(phi0_2) ;
239
240 // Vector that goes from center1 to center2
241 const double deltax = x0_2 - x0_1 ;
242 const double deltay = y0_2 - y0_1 ;
243
244 // Intersections of the circles, can be at most two
245 double phi1[2] ;
246 double phi2[2] ;
247 int nsolutions = 1;
248
249 // The phi of the delta vector.
250 const double phi = - atan2(deltax, deltay) ;
251 const double phinot = phi > 0 ? phi - TMath::Pi() : phi + TMath::Pi() ;
252 phi1[0] = r_1 < 0 ? phi : phinot ;
253 phi2[0] = r_2 > 0 ? phi : phinot ;
254
255 // These radii do NOT have a sign instead
256 const double R1 = fabs(r_1) ;
257 const double R2 = fabs(r_2) ;
258 const double Rmin = R1 < R2 ? R1 : R2 ;
259 const double Rmax = R1 > R2 ? R1 : R2 ;
260 const double dX = hypot(deltax, deltay) ;
261
262 if (!parallel && dX + Rmin > Rmax && dX < R1 + R2) {
263 // Circles intersect in two points
264 nsolutions = 2 ;
265
266 // This is just the law of cosines
267 const double ddphi1 = acos((dX * dX - R2 * R2 + R1 * R1) / (2.*dX * R1)) ;
268 phi1[1] = phidomain(phi1[0] + ddphi1) ;
269 phi1[0] = phidomain(phi1[0] - ddphi1) ;
270
271 const double ddphi2 = acos((dX * dX - R1 * R1 + R2 * R2) / (2.*dX * R2)) ;
272 phi2[1] = phidomain(phi2[0] - ddphi2) ;
273 phi2[0] = phidomain(phi2[0] + ddphi2) ;
274
275 } else if (dX < Rmax) {
276 // Tangent or non-intersecting circles, one inside the other (only one POCA)
277 if (R1 > R2) phi2[0] = r_2 < 0 ? phi : phinot ;
278 else phi1[0] = r_1 < 0 ? phi : phinot ;
279 }
280 // else: tangent or non-intersecting circles, outside of each other (only one POCA)
281 // what we saved in phi1 and phi2 gives already the correct solution
282
283 // Intersections of the circles (cartesian)
284 double x1[2], y1[2], x2[2], y2[2];
285 for (int i = 0; i < nsolutions; i++) {
286 x1[i] = r_1 * sin(phi1[i]) + x0_1 ;
287 y1[i] = -r_1 * cos(phi1[i]) + y0_1 ;
288 x2[i] = r_2 * sin(phi2[i]) + x0_2 ;
289 y2[i] = -r_2 * cos(phi2[i]) + y0_2 ;
290 }
291
292 // 2) Find the best solution for z by running multiples of 2pi from the xy intersection(s)
293 double z1(0), z2(0);
294 bool first = true;
295 int ibest = 0;
296 const int nturnsmax = 10; // Max number of turns we try backwards and forwards
297
298 // Loop on all xy-plane solutions
299 for (int i = 0; i < nsolutions; ++i) {
300 const double l1 = helix1.getArcLength2DAtXY(x1[i], y1[i]);
301 const double l2 = helix2.getArcLength2DAtXY(x2[i], y2[i]);
302
303 // Loop on helix1 turns, save corresponding z positions
304 std::vector<double> z1s;
305 for (int n1 = 0; n1 <= nturnsmax; ++n1) {
306 bool added = false;
307 // Try forwards and backwards
308 for (int sn1 : {n1, -n1}) {
309 const double tmpz1 = helix1.getPositionAtArcLength2D(l1 + sn1 * TMath::TwoPi() / omega_1).Z();
310 if (sn1 == 0 || (-82 <= tmpz1 && tmpz1 <= 158)) {
311 // Only keep the 0th turn and those inside CDC volume
312 z1s.push_back(tmpz1);
313 added = true;
314 }
315 if (sn1 == 0)
316 break; // Do not store 0th turn twice
317 }
318 // If we did not add any point we are already outside CDC volume both backwards and forwards
319 if (!added)
320 break;
321 }
322
323 // Loop on helix2 turns, find closest approach to one of helix1 points
324 for (int n2 = 0; n2 <= nturnsmax; ++n2) {
325 bool tried = false;
326 // Try forwards and backwards
327 for (int sn2 : {n2, -n2}) {
328 const double tmpz2 = helix2.getPositionAtArcLength2D(l2 + sn2 * TMath::TwoPi() / omega_2).Z();
329 if (sn2 == 0 || (-82 <= tmpz2 && tmpz2 <= 158)) {
330 // Only keep the 0th turn and those inside CDC volume
331 tried = true;
332 // Find the tmpz1 closest to tmpz2
333 const auto i1best = std::min_element(
334 z1s.cbegin(), z1s.cend(), [&tmpz2](const double & z1a, const double & z1b) {
335 return fabs(z1a - tmpz2) < fabs(z1b - tmpz2);
336 });
337 const double tmpz1 = *i1best;
338 // Keep the solution where the z distance of closest approach is minimum
339 if (first || fabs(tmpz1 - tmpz2) < fabs(z1 - z2)) {
340 ibest = i;
341 first = false;
342 z1 = tmpz1;
343 z2 = tmpz2;
344 flt1 = l1;
345 flt2 = l2;
346 }
347 }
348 if (n2 == 0)
349 break; // Do not try 0th turn twice
350 }
351 // If we did not try any point we are already outside CDC volume both backwards and forwards
352 if (!tried)
353 break;
354 }
355 }
356
357 vertex.SetX(0.5 * (x1[ibest] + x2[ibest]));
358 vertex.SetY(0.5 * (y1[ibest] + y2[ibest]));
359 vertex.SetZ(0.5 * (z1 + z2));
360
361 return hypot(x2[ibest] - x1[ibest], y2[ibest] - y1[ibest], z2 - z1);
362 }
ROOT::Math::XYZVector getPositionAtArcLength2D(const double &arcLength2D) const
Calculates the position on the helix at the given two dimensional arc length.
Definition: Helix.cc:201

◆ phidomain()

double phidomain ( const double  phi)
static

the domain of phi

Definition at line 205 of file HelixUtils.cc.

206 {
207 double rc = phi ;
208 if (phi < -TMath::Pi()) rc += TMath::TwoPi();
209 else if (phi > TMath::Pi()) rc -= TMath::TwoPi();
210 return rc ;
211 }

◆ printVertexPar()

void printVertexPar ( const Belle2::B2Vector3D position,
const Belle2::B2Vector3D momentum,
int  charge 
)
static

Print the vertex parameters.

Definition at line 120 of file HelixUtils.cc.

121 {
122 for (int i = 0; i < 3; ++i)
123 B2INFO(vertexParName(i + 1).c_str() << position[i]);
124 for (int i = 0; i < 3; ++i)
125 B2INFO(vertexParName(i + 4).c_str() << momentum[i]);
126 B2INFO("charge: " << charge);
127
128 }
static std::string vertexParName(int i)
map of the vertex parameters by list index
Definition: HelixUtils.cc:106

◆ vertexFromHelix()

void vertexFromHelix ( const Belle2::Helix helix,
double  L,
double  Bz,
ROOT::Math::XYZVector &  position,
ROOT::Math::XYZVector &  momentum,
int &  charge 
)
static

helix --> vertex

Definition at line 23 of file HelixUtils.cc.

27 {
28 position = helix.getPositionAtArcLength2D(L);
29 momentum = helix.getMomentumAtArcLength2D(L, Bz);
30 charge = helix.getChargeSign();
31 }
short getChargeSign() const
Return track charge sign (1, 0 or -1).
Definition: Helix.cc:114
ROOT::Math::XYZVector getMomentumAtArcLength2D(const double &arcLength2D, const double &bz) const
Calculates the momentum vector at the given two dimensional arc length.
Definition: Helix.cc:270

◆ vertexParName()

std::string vertexParName ( int  i)
static

map of the vertex parameters by list index

Definition at line 106 of file HelixUtils.cc.

107 {
108 std::string rc ;
109 switch (i) {
110 case 1 : rc = "x : " ; break ;
111 case 2 : rc = "y : " ; break ;
112 case 3 : rc = "z : " ; break ;
113 case 4 : rc = "px : " ; break ;
114 case 5 : rc = "py : " ; break ;
115 case 6 : rc = "pz : " ; break ;
116 }
117 return rc ;
118 }

The documentation for this class was generated from the following files: