Skip to content

Commit

Permalink
Initialize_Canonical and getDistanceToHelix updates for reference poi…
Browse files Browse the repository at this point in the history 
…nts other than (0,0,0)
  • Loading branch information
@tmadlener
Jan Klamka authored and tmadlener committed Oct 19, 2022
1 parent eda5bc0 commit ce93d98
Show file tree
Hide file tree
Showing 2 changed files with 136 additions and 32 deletions.
28 changes: 27 additions & 1 deletion source/include/HelixClassT.h
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
#ifndef HELIXCLASST_H
#define HELIXCLASST_H 1


#include <vector>
#include <cmath>

Expand Down Expand Up @@ -105,6 +104,23 @@ class HelixClassT {
*/
void Initialize_Canonical(FloatT phi0, FloatT d0, FloatT z0, FloatT omega,
FloatT tanlambda, FloatT B);

/**
* Canonical (LEP-wise) parameterisation with the following parameters<br>
* - phi0 - Phi angle of momentum vector at the point of<br>
* closest approach to IP in R-Phi plane;
* - d0 - signed distance of closest approach in R-Phi plane;<br>
* - z0 - z coordinate of the point of closest approach to IP
* in R-Phi plane;<br>
* - omega - signed curvature;<br>
* - tanlambda - tangent of dip angle;<br>
* - B - magnetic field (in Tesla)<br>
* - pos - reference point at which the helix parameteres are given
* (by default PCA is used as the reference point);<br>
*/
void Initialize_Canonical(FloatT phi0, FloatT d0, FloatT z0,
FloatT omega, FloatT tanLambda, FloatT B, FloatT * pos);

/**
* Returns momentum of particle at the point of closest approach <br>
* to IP <br>
Expand Down Expand Up @@ -160,6 +176,16 @@ class HelixClassT {
*/
FloatT getYC() const { return _yCentre; }

/**
* Returns x coordinate of point of closest approach
*/
FloatT getXPCA() const { return _xAtPCA; }

/**
* Returns y coordinate of point of closest approach
*/
FloatT getYPCA() const { return _yAtPCA; }


/**
* Returns radius of circumference
Expand Down
140 changes: 109 additions & 31 deletions source/include/HelixClassT.ipp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,8 @@
#include <iostream>
#include "ced_cli.h"

#include <streamlog/streamlog.h>

template<typename FloatT>
void HelixClassT<FloatT>::Initialize_VP(FloatT * pos, FloatT * mom, FloatT q, FloatT B) {
_referencePoint[0] = pos[0];
Expand Down Expand Up @@ -112,6 +114,36 @@ void HelixClassT<FloatT>::Initialize_Canonical(FloatT phi0, FloatT d0, FloatT z0
_bField = B;
}

template<typename FloatT>
void HelixClassT<FloatT>::Initialize_Canonical(FloatT phi0, FloatT d0, FloatT z0,
FloatT omega, FloatT tanLambda, FloatT B, FloatT* pos) {
_omega = omega;
_d0 = d0;
_phi0 = phi0;
_z0 = z0;
_tanLambda = tanLambda;
_charge = omega/fabs(omega);
_radius = 1./fabs(omega);
_referencePoint[0] = pos[0];
_referencePoint[1] = pos[1];
_referencePoint[2] = pos[2];
_xAtPCA = -_d0*sin(_phi0) + _referencePoint[0];
_yAtPCA = _d0*cos(_phi0) + _referencePoint[1];
_pxy = _FCT*B*_radius;
_momentum[0] = _pxy*cos(_phi0);
_momentum[1] = _pxy*sin(_phi0);
_momentum[2] = _tanLambda * _pxy;
_pxAtPCA = _momentum[0];
_pyAtPCA = _momentum[1];
_phiMomRefPoint = atan2(_momentum[1],_momentum[0]);
_xCentre = _referencePoint[0] +
_radius*cos(_phi0-_const_pi2*_charge);
_yCentre = _referencePoint[1] +
_radius*sin(_phi0-_const_pi2*_charge);
_phiAtPCA = atan2(-_yCentre,-_xCentre);
_phiRefPoint = _phiAtPCA ;
_bField = B;
}

template<typename FloatT>
void HelixClassT<FloatT>::Initialize_BZ(FloatT xCentre, FloatT yCentre, FloatT radius,
Expand Down Expand Up @@ -248,10 +280,9 @@ FloatT HelixClassT<FloatT>::getPointInXY(FloatT x0, FloatT y0, FloatT ax, FloatT
tt2 = -_charge*dphi2*_radius/_pxy;

if (tt1 < 0. )
std::cout << "WARNING " << tt1 << std::endl;
streamlog_out(WARNING) << "Time : " << tt1 << " < 0!" << std::endl;
if (tt2 < 0. )
std::cout << "WARNING " << tt2 << std::endl;

streamlog_out(WARNING) << "Time : " << tt2 << " < 0!" << std::endl;

if (tt1 < tt2) {
point[0] = xx1;
Expand Down Expand Up @@ -341,9 +372,9 @@ FloatT HelixClassT<FloatT>::getPointOnCircle(FloatT Radius, FloatT * ref, FloatT
tt2 = -_charge*dphi2*_radius/_pxy;

if (tt1 < 0. )
std::cout << "WARNING " << tt1 << std::endl;
streamlog_out(WARNING) << "Time : " << tt1 << " < 0!" << std::endl;
if (tt2 < 0. )
std::cout << "WARNING " << tt2 << std::endl;
streamlog_out(WARNING) << "Time : " << tt2 << " < 0!" << std::endl;


FloatT time2;
Expand Down Expand Up @@ -594,9 +625,6 @@ FloatT HelixClassT<FloatT>::getDistanceToHelix(HelixClassT * helix, FloatT * pos
FloatT xSect2 = x02 + rad2*cos(phi2);
FloatT ySect2 = y02 + rad2*sin(phi2);

// std::cout << "(xSect1,ySect1)=(" << xSect1 << "," << ySect1 << ")" << std::endl;
// std::cout << "(xSect2,ySect2)=(" << xSect2 << "," << ySect2 << ")" << std::endl;

FloatT temp21[3];
FloatT temp22[3];

Expand All @@ -608,25 +636,52 @@ FloatT HelixClassT<FloatT>::getDistanceToHelix(HelixClassT * helix, FloatT * pos
FloatT charge2 = helix->getCharge();
FloatT pxy2 = helix->getPXY();
FloatT pz2 = helix->getMomentum()[2];
if (deltaPhi21 < 0 && charge2 < 0) {
deltaPhi21 += _const_2pi;
}
else if (deltaPhi21 > 0 && charge2 > 0) {
deltaPhi21 -= _const_2pi;
}

if (deltaPhi22 < 0 && charge2 < 0) {
deltaPhi22 += _const_2pi;
}
else if (deltaPhi22 > 0 && charge2 > 0) {
deltaPhi22 -= _const_2pi;
float xAtPCA2 = helix->getXPCA();
float yAtPCA2 = helix->getYPCA();

// if helix initialized with no input ref. point or with (0,0,0)
if ( (ref2[0]==xAtPCA2 && ref2[1]==yAtPCA2 && ref2[2]==helix->getZ0())
|| (ref2[0]==0 && ref2[1]==0 && ref2[2]==0) ) {
if (deltaPhi21 < 0 && charge2 < 0) {
deltaPhi21 += _const_2pi;
}
else if (deltaPhi21 > 0 && charge2 > 0) {
deltaPhi21 -= _const_2pi;
}

if (deltaPhi22 < 0 && charge2 < 0) {
deltaPhi22 += _const_2pi;
}
else if (deltaPhi22 > 0 && charge2 > 0) {
deltaPhi22 -= _const_2pi;
}
} else { // any other ref. point
if (fabs(deltaPhi21) > M_PI) {
if (deltaPhi21 < 0 && charge2 > 0) {
deltaPhi21 += _const_2pi;
}
else if (deltaPhi21 > 0 && charge2 < 0) {
deltaPhi21 -= _const_2pi;
}
}
if (fabs(deltaPhi22) > M_PI) {
if (deltaPhi22 < 0 && charge2 > 0) {
deltaPhi22 += _const_2pi;
}
else if (deltaPhi22 > 0 && charge2 < 0) {
deltaPhi22 -= _const_2pi;
}
}
}

FloatT time21 = -charge2*deltaPhi21*rad2/ pxy2;
FloatT time22 = -charge2 * deltaPhi22 * rad2 / pxy2;

FloatT Z21 = ref2[2] + time21 * pz2;
FloatT Z22 = ref2[2] + time22 * pz2;
FloatT tanLambda2 = helix->getTanLambda();

FloatT Z21 = ref2[2] - charge2*deltaPhi21*rad2*tanLambda2;
FloatT Z22 = ref2[2] - charge2*deltaPhi22*rad2*tanLambda2;

temp21[0] = xSect1;
temp21[1] = ySect1;
Expand All @@ -650,10 +705,35 @@ FloatT HelixClassT<FloatT>::getDistanceToHelix(HelixClassT * helix, FloatT * pos
FloatT charge1 = _charge;
FloatT pxy1 = _pxy;
FloatT pz1 = _momentum[2];
if (deltaPhi11 < 0 && charge1 < 0) {
deltaPhi11 += _const_2pi;
} else if (deltaPhi11 > 0 && charge1 > 0) {
deltaPhi11 -= _const_2pi;

FloatT xAtPCA1 = _xAtPCA;
FloatT yAtPCA1 = _yAtPCA;

// if helix initialized with no input ref. point or with (0,0,0)
if ( (ref1[0]==xAtPCA1 && ref1[1]==yAtPCA1 && ref1[2]==_z0)
|| (ref1[0]==0 && ref1[1]==0 && ref1[2]==0) ) {
if (deltaPhi11 < 0 && charge1 < 0) {
deltaPhi11 += _const_2pi;
} else if (deltaPhi11 > 0 && charge1 > 0) {
deltaPhi11 -= _const_2pi;
}
} else { // any other ref. point
if (fabs(deltaPhi11) > M_PI) {
if (deltaPhi11 < 0 && charge1 > 0) {
deltaPhi11 += _const_2pi;
}
else if (deltaPhi11 > 0 && charge1 < 0) {
deltaPhi11 -= _const_2pi;
}
}
if (fabs(deltaPhi12) > M_PI) {
if (deltaPhi12 < 0 && charge1 > 0) {
deltaPhi12 += _const_2pi;
}
else if (deltaPhi12 > 0 && charge1 < 0) {
deltaPhi12 -= _const_2pi;
}
}
}

if (deltaPhi12 < 0 && charge1 < 0) {
Expand All @@ -665,8 +745,10 @@ FloatT HelixClassT<FloatT>::getDistanceToHelix(HelixClassT * helix, FloatT * pos
FloatT time11 = -charge1 * deltaPhi11 * rad1 / pxy1;
FloatT time12 = -charge1 * deltaPhi12 * rad1 / pxy1;

FloatT Z11 = ref1[2] + time11 * pz1;
FloatT Z12 = ref1[2] + time12 * pz1;
FloatT tanLambda1 = getTanLambda();

FloatT Z11 = ref1[2] - charge1*deltaPhi11*rad1*tanLambda1;
FloatT Z12 = ref1[2] - charge1*deltaPhi12*rad1*tanLambda1;

temp11[0] = xSect1;
temp11[1] = ySect1;
Expand All @@ -675,10 +757,6 @@ FloatT HelixClassT<FloatT>::getDistanceToHelix(HelixClassT * helix, FloatT * pos
temp12[1] = ySect2;
temp12[2] = Z12;

// std::cout << "temp11 = " << temp11[0] << " " << temp11[1] << " " <<
// temp11[2] << std::endl; std::cout << "temp12 = " << temp12[0] << " "
// << temp12[1] << " " << temp12[2] << std::endl;

FloatT Dist1 = 0;
FloatT Dist2 = 0;

Expand Down

0 comments on commit ce93d98

Please sign in to comment.