From 62f012ffe60dee3b26cce38dcf904d7bdaa5ba52 Mon Sep 17 00:00:00 2001 From: Jan Klamka Date: Mon, 19 Sep 2022 20:19:03 +0200 Subject: [PATCH] Add getXYPCA() and description of getDistanceToHelix(), remove commented printout --- source/include/HelixClass.h | 31 +++++++++++++--- source/src/HelixClass.cc | 73 +++++++++++-------------------------- 2 files changed, 46 insertions(+), 58 deletions(-) diff --git a/source/include/HelixClass.h b/source/include/HelixClass.h index 73ad858..1ddb73f 100644 --- a/source/include/HelixClass.h +++ b/source/include/HelixClass.h @@ -113,9 +113,9 @@ class HelixClass { * in R-Phi plane;
* - omega - signed curvature;
* - tanlambda - tangent of dip angle;
- * - B - magnetic field (in Tesla)
- * - pos - reference point at which the helix parameteres are given - * (by default PCA is used as the reference point);
+ * - B - magnetic field (in Tesla);
+ * - pos - reference point at which the helix parameteres are given + * (by default PCA is used as the reference point)
*/ void Initialize_Canonical(float phi0, float d0, float z0, float omega, float tanLambda, float B, float * pos); @@ -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 @@ -242,10 +251,20 @@ class HelixClass { */ float getPointOnCircle(float Radius, float * ref, float * point); - /** Returns distance between two helixes
+ /** + * 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.
+ * Input :
+ * helix - the method should be called on the helix with higher transverse + * momentum, with the other helix provided as an input
* Output :
- * pos[3] - position of the point of closest approach
- * mom[3] - momentum of V0
+ * pos[3] - position of the point exactly between the points of closest + * approach (vertex)
+ * mom[3] - sum of the track momenta at the vertex
*/ float getDistanceToHelix(HelixClass * helix, float * pos, float * mom); diff --git a/source/src/HelixClass.cc b/source/src/HelixClass.cc index 85daa2d..80a96a6 100644 --- a/source/src/HelixClass.cc +++ b/source/src/HelixClass.cc @@ -1,9 +1,13 @@ #include "HelixClass.h" + #include #include #include #include "ced_cli.h" +// ----- include for verbosity dependend logging --------- +#include "marlin/VerbosityLevels.h" + HelixClass::HelixClass() {} HelixClass::~HelixClass() {} @@ -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; @@ -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; } @@ -235,6 +229,14 @@ float HelixClass::getYC() { return _yCentre; } +float HelixClass::getXPCA() { + return _xAtPCA; +} + +float HelixClass::getYPCA() { + return _yAtPCA; +} + float HelixClass::getRadius() { return _radius; } @@ -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) { @@ -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; @@ -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; @@ -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]; @@ -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()) @@ -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; @@ -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]; @@ -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) @@ -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; @@ -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