Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Initialize_Canonical() and getDistanceToHelix() updates for reference points other than (0,0,0) #29

Open
wants to merge 7 commits into
base: master
Choose a base branch
from
Open

Initialize_Canonical() and getDistanceToHelix() updates for reference points other than (0,0,0) #29

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
ce93d98
Initialize_Canonical and getDistanceToHelix updates for reference poi…
Sep 6, 2022
1ed3f93
Minor fixes after rebasing
Oct 23, 2022
5d931ee
Merge branch 'iLCSoft:master' into HelixClassFix
jfklama Oct 23, 2022
633ad6f
Merge branch 'HelixClassFix' of https://github.com/jfklama/MarlinUtil…
Oct 23, 2022
75a0e92
Unit tests for getDistanceToHelix() and minor change to the function
Nov 9, 2022
c817eb9
Leave only unit tests relevant to the most recent HelixClass changes
Nov 15, 2022
179af9c
Fix _xCentre, _yCentre in the new Initialize_Canonical
Nov 16, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 27 additions & 1 deletion source/include/HelixClassT.h
View file
Open in desktop
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 @@ -108,6 +107,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 @@ -163,6 +179,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
138 changes: 101 additions & 37 deletions source/include/HelixClassT.ipp
View file
Open in desktop
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(_referencePoint[1]-_yCentre,_referencePoint[0]-_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,46 @@ 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;
float xAtPCA2 = helix->getXPCA();
float yAtPCA2 = helix->getYPCA();


// if helix initialized with no input ref. point, for backward compatibility
if ( ref2[0]==xAtPCA2 && ref2[1]==yAtPCA2 && ref2[2]==helix->getZ0() ) {
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 if (deltaPhi22 > 0 && charge2 > 0) {
deltaPhi22 -= _const_2pi;

else { // for non-zero ref. points given as an input
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 @@ -635,9 +684,6 @@ FloatT HelixClassT<FloatT>::getDistanceToHelix(HelixClassT * helix, FloatT * pos
temp22[1] = ySect2;
temp22[2] = Z22;

// std::cout << "temp21 = " << temp21[0] << " " << temp21[1] << " " <<
// temp21[2] << std::endl; std::cout << "temp22 = " << temp22[0] << " "
// << temp22[1] << " " << temp22[2] << std::endl;

FloatT temp11[3];
FloatT temp12[3];
Expand All @@ -650,23 +696,45 @@ 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, for backward compatibility
if ( ref1[0]==xAtPCA1 && ref1[1]==yAtPCA1 && ref1[2]==_z0 ) {
if (deltaPhi11 < 0 && charge1 < 0)
deltaPhi11 += _const_2pi;
else if (deltaPhi11 > 0 && charge1 > 0)
deltaPhi11 -= _const_2pi;
if (deltaPhi12 < 0 && charge1 < 0)
deltaPhi12 += _const_2pi;
else if (deltaPhi12 > 0 && charge1 > 0)
deltaPhi12 -= _const_2pi;
}

if (deltaPhi12 < 0 && charge1 < 0) {
deltaPhi12 += _const_2pi;
} else if (deltaPhi12 > 0 && charge1 > 0) {
deltaPhi12 -= _const_2pi;
else { // for non-zero ref. points given as an input
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;
}
}


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 +743,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
Loading