Skip to content

Commit

Permalink
Add getXYPCA() and description of getDistanceToHelix(), remove commen…
Browse files Browse the repository at this point in the history 
…ted printout
  • Loading branch information
Jan Klamka committed Sep 19, 2022
1 parent 05ff212 commit 62f012f
Show file tree
Hide file tree
Showing 2 changed files with 46 additions and 58 deletions.
31 changes: 25 additions & 6 deletions source/include/HelixClass.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,9 +113,9 @@ class HelixClass {
* 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>
* - 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(float phi0, float d0, float z0,
float omega, float tanLambda, float B, float * pos);
Expand Down Expand Up @@ -175,6 +175,15 @@ class HelixClass {
*/
float getYC();

/**
* Returns x coordinate of point of closest approach
*/
float getXPCA();

/**
* Returns y coordinate of point of closest approach
*/
float getYPCA();

/**
* Returns radius of circumference
Expand Down Expand Up @@ -242,10 +251,20 @@ class HelixClass {
*/
float getPointOnCircle(float Radius, float * ref, float * point);

/** Returns distance between two helixes <br>
/**
* Returns distance between two helixes. The two closest
* points on the helices are found using projections to the xy plane and
* getPointOnCircle() method. If they intersect (more than once), two
* candidates for the intersection points are calculated from the parameters
* of each helix. The smaller distance between the corresponding points on
* the two helices is returned. <br>
* Input : <br>
* helix - the method should be called on the helix with higher transverse
* momentum, with the other helix provided as an input <br>
* Output : <br>
* pos[3] - position of the point of closest approach <br>
* mom[3] - momentum of V0 <br>
* pos[3] - position of the point exactly between the points of closest
* approach (vertex) <br>
* mom[3] - sum of the track momenta at the vertex <br>
*/
float getDistanceToHelix(HelixClass * helix, float * pos, float * mom);

Expand Down
73 changes: 21 additions & 52 deletions source/src/HelixClass.cc
View file
Original file line number Diff line number Diff line change
@@ -1,9 +1,13 @@
#include "HelixClass.h"

#include <math.h>
#include <stdlib.h>
#include <iostream>
#include "ced_cli.h"

// ----- include for verbosity dependend logging ---------
#include "marlin/VerbosityLevels.h"

HelixClass::HelixClass() {}

HelixClass::~HelixClass() {}
Expand Down Expand Up @@ -48,16 +52,6 @@ void HelixClass::Initialize_VP(float * pos, float * mom, float q, float B) {

_d0 = float(d0);

// if (fabs(_d0)>0.001 ) {
// std::cout << "New helix : " << std::endl;
// std::cout << " Position : " << pos[0]
// << " " << pos[1]
// << " " << pos[2] << std::endl;
// std::cout << " Radius = " << _radius << std::endl;
// std::cout << " RC = " << sqrt(_xCentre*_xCentre+_yCentre*_yCentre) << std::endl;
// std::cout << " D0 = " << _d0 << std::endl;
// }

_pxAtPCA = _pxy*cos(_phi0);
_pyAtPCA = _pxy*sin(_phi0);
float deltaPhi = _phiRefPoint - _phiAtPCA;
Expand Down Expand Up @@ -140,7 +134,7 @@ void HelixClass::Initialize_Canonical(float phi0, float d0, float z0,
_radius*cos(_phi0-_const_pi2*_charge);
_yCentre = _referencePoint[1] +
_radius*sin(_phi0-_const_pi2*_charge);
_phiAtPCA = atan2(-_yCentre,-_xCentre);
_phiAtPCA = atan2(_referencePoint[1]-_yCentre,_referencePoint[0]-_xCentre);
_phiRefPoint = _phiAtPCA ;
_bField = B;
}
Expand Down Expand Up @@ -235,6 +229,14 @@ float HelixClass::getYC() {
return _yCentre;
}

float HelixClass::getXPCA() {
return _xAtPCA;
}

float HelixClass::getYPCA() {
return _yAtPCA;
}

float HelixClass::getRadius() {
return _radius;
}
Expand Down Expand Up @@ -322,9 +324,9 @@ float HelixClass::getPointInXY(float x0, float y0, float ax, float ay,
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) {
Expand Down Expand Up @@ -414,9 +416,9 @@ float HelixClass::getPointOnCircle(float Radius, float * ref, float * point) {
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;


float time2;
Expand Down Expand Up @@ -598,18 +600,9 @@ float HelixClass::getDistanceToHelix(HelixClass * helix, float * pos, float * mo
float x02 = helix->getXC();
float y02 = helix->getYC();

//std::cout << "rad1: " << getRadius() << ", phi01: " << _phi0 << ", charge: " << _charge << ", sin(phi0-M_PI*charge): " << sin(_phi0-M_PI*_charge) << ", pos[1] + sin(phi0-M_PI*charge): " << getReferencePoint()[1] + sin(_phi0-M_PI*_charge) << std::endl;
//std::cout << "rad2: " << helix->getRadius() << ", phi01: " << helix->getPhi0() << ", charge: " << helix->getCharge() << ", sin(phi0-M_PI*charge): " << sin(helix->getPhi0()-M_PI*helix->getCharge()) << ", pos[1] + sin(phi0-M_PI*charge): " << helix->getReferencePoint()[1] + sin(helix->getPhi0()-M_PI*helix->getCharge()) << std::endl;

//std::cout << "(x01,y01)=(" << x01 << "," << y01 << ")" << std::endl;
//std::cout << "(x02,y02)=(" << x02 << "," << y02 << ")" << std::endl;

float rad1 = getRadius();
float rad2 = helix->getRadius();

//std::cout << "rad1=" << rad1 << std::endl;
//std::cout << "rad2=" << rad2 << std::endl;

float distance = sqrt((x01-x02)*(x01-x02)+(y01-y02)*(y01-y02));

bool singlePoint = true;
Expand Down Expand Up @@ -672,9 +665,6 @@ float HelixClass::getDistanceToHelix(HelixClass * helix, float * pos, float * mo
float xSect2 = x02 + rad2*cos(phi2);
float ySect2 = y02 + rad2*sin(phi2);

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

float temp21[3];
float temp22[3];

Expand All @@ -687,8 +677,8 @@ float HelixClass::getDistanceToHelix(HelixClass * helix, float * pos, float * mo
float pxy2 = helix->getPXY();
float pz2 = helix->getMomentum()[2];

float xAtPCA2 = -helix->getD0()*sin(helix->getPhi0());
float yAtPCA2 = helix->getD0()*cos(helix->getPhi0());
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())
Expand Down Expand Up @@ -725,14 +715,10 @@ float HelixClass::getDistanceToHelix(HelixClass * helix, float * pos, float * mo
}
}
}
//std::cout << "phiRef2: " << phiI2 << ", phiSect21: " << phiF21 << ", phiSect22: " << phiF22 << ", deltaPhi21: " << deltaPhi21 << ", deltaPhi22: " << deltaPhi22 << std::endl;

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

//float Z21 = ref2[2] + time21*pz2;
//float Z22 = ref2[2] + time22*pz2;

float tanLambda2 = helix->getTanLambda();

float Z21 = ref2[2] - charge2*deltaPhi21*rad2*tanLambda2;
Expand All @@ -742,10 +728,6 @@ float HelixClass::getDistanceToHelix(HelixClass * helix, float * pos, float * mo
temp22[0] = xSect2; 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;


float temp11[3];
float temp12[3];

Expand All @@ -758,8 +740,8 @@ float HelixClass::getDistanceToHelix(HelixClass * helix, float * pos, float * mo
float pxy1 = _pxy;
float pz1 = _momentum[2];

float xAtPCA1 = -_d0*sin(_phi0);
float yAtPCA1 = _d0*cos(_phi0);
float xAtPCA1 = _xAtPCA;
float 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)
Expand Down Expand Up @@ -797,28 +779,18 @@ float HelixClass::getDistanceToHelix(HelixClass * helix, float * pos, float * mo
}
}
}
//std::cout << "phiRef1: " << phiI1 << ", phiSect11: " << phiF11 << ", phiSect12: " << phiF12 << ", deltaPhi11: " << deltaPhi11 << ", deltaPhi12: " << deltaPhi12 << std::endl;

float time11 = -charge1*deltaPhi11*rad1/pxy1;
float time12 = -charge1*deltaPhi12*rad1/pxy1;

//float Z11 = ref1[2] + time11*pz1;
//float Z12 = ref1[2] + time12*pz1;

float tanLambda1 = getTanLambda();

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

//std::cout << "rad1: " << getRadius() << ", deltaPhi11: " << deltaPhi11 << ", charge1: " << charge1 << ", tanLambda1: " << tanLambda1 << ", charge1*deltaPhi11*rad1*tanLambda1: " << charge1*deltaPhi11*rad1*tanLambda1 << ", pos1[2] - charge1*deltaPhi11*rad1*tanLambda1: " << getReferencePoint()[2] - charge1*deltaPhi11*rad1*tanLambda1 << std::endl;
//std::cout << "rad2: " << helix->getRadius() << ", deltaPhi21: " << deltaPhi21 << ", charge2: " << charge2 << ", tanLambda2: " << tanLambda2 << ", charge2*deltaPhi21*rad2*tanLambda2: " << charge2*deltaPhi21*rad2*tanLambda2 << ", pos2[2] - charge2*deltaPhi21*rad2*tanLambda2: " << helix->getReferencePoint()[2] - charge2*deltaPhi21*rad2*tanLambda2 << std::endl;

temp11[0] = xSect1; temp11[1] = ySect1; temp11[2] = Z11;
temp12[0] = xSect2; 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;

float Dist1 = 0;
float Dist2 = 0;

Expand All @@ -827,9 +799,6 @@ float HelixClass::getDistanceToHelix(HelixClass * helix, float * pos, float * mo
Dist2 += (temp12[j]-temp22[j])*(temp12[j]-temp22[j]);
}

//std::cout << "Dist1 = " << Dist1 << std::endl;
//std::cout << "Dist2 = " << Dist2 << std::endl;

if (Dist1<Dist2) {
for (int l=0;l<3;++l) {
pos1[l] = temp11[l];
Expand Down

0 comments on commit 62f012f

Please sign in to comment.