From f17db08fc8911c6da9b7927f7e560bdbaed98114 Mon Sep 17 00:00:00 2001 From: Xin Dong Date: Wed, 1 Oct 2025 19:11:46 -0400 Subject: [PATCH 01/17] added SecondaryVertex factory using Helix functions --- src/algorithms/reco/Helix.cc | 659 ++++++++++++++++++ src/algorithms/reco/Helix.h | 238 +++++++ src/algorithms/reco/SecondaryVertices.cc | 125 ++++ src/algorithms/reco/SecondaryVertices.h | 37 + src/algorithms/reco/SecondaryVerticesConfig.h | 22 + .../reco/SecondaryVertices_factory.h | 50 ++ src/global/reco/reco.cc | 6 +- src/services/io/podio/JEventProcessorPODIO.cc | 1 + 8 files changed, 1137 insertions(+), 1 deletion(-) create mode 100644 src/algorithms/reco/Helix.cc create mode 100644 src/algorithms/reco/Helix.h create mode 100644 src/algorithms/reco/SecondaryVertices.cc create mode 100644 src/algorithms/reco/SecondaryVertices.h create mode 100644 src/algorithms/reco/SecondaryVerticesConfig.h create mode 100644 src/factories/reco/SecondaryVertices_factory.h diff --git a/src/algorithms/reco/Helix.cc b/src/algorithms/reco/Helix.cc new file mode 100644 index 0000000000..eed158857e --- /dev/null +++ b/src/algorithms/reco/Helix.cc @@ -0,0 +1,659 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Xin Dong, Rongrong Ma + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include "algorithms/reco/Helix.h" +#include +#include +#include + +namespace eicrecon { + +const double Helix::NoSolution = 3.e+33; + +// basic constructor +Helix::Helix(double c, double d, double phase, const edm4hep::Vector3f& o, int h) +{ + setParameters(c, d, phase, o, h); +} + +// momentum in GeV, position in cm, B in Tesla +Helix::Helix(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q) +{ + setParameters(p, o, B, q); +} + +// construct using ReconstructParticle +Helix::Helix(const edm4eic::ReconstructedParticle& p, const double b_field) { + const auto& tracks = p.getTracks(); + for (const auto& trk : tracks) { + const auto& traj = trk.getTrajectory(); + const auto& trkPars = traj.getTrackParameters(); + for (const auto& par : trkPars) { + setParameters(par, b_field); + } + } +} + +// construct using TrackParameteters +void Helix::setParameters(const edm4eic::TrackParameters& trk, const double b_field) { + const auto mom = edm4hep::utils::sphericalToVector(1.0 / std::abs(trk.getQOverP()), trk.getTheta(), trk.getPhi()); + const auto charge = std::copysign(1., trk.getQOverP()); + const auto phi = trk.getPhi(); + const auto loc = trk.getLoc(); + edm4hep::Vector3f pos( -1. * loc.a * sin(phi), loc.a * cos(phi), loc.b); // PCA point + + setParameters(mom*edm4eic::unit::GeV/dd4hep::GeV, pos*edm4eic::unit::mm/edm4eic::unit::cm, b_field, charge); +} + +void Helix::setParameters(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q) +{ + mH = (q*B <= 0) ? 1 : -1; + if(p.y == 0 && p.x == 0) + setPhase((M_PI/4)*(1-2.*mH)); + else + setPhase(atan2(p.y,p.x)-mH*M_PI/2); + setDipAngle(atan2(p.z,edm4hep::utils::magnitudeTransverse(p))); + mOrigin = o; + + double curvature_val = fabs((dd4hep::c_light*dd4hep::nanosecond/dd4hep::meter*q*B/dd4hep::tesla)/ + (edm4hep::utils::magnitude(p)/dd4hep::GeV*mCosDipAngle)/dd4hep::meter); + setCurvature(curvature_val); +} + +void Helix::setParameters(double c, double dip, double phase, + const edm4hep::Vector3f& o, int h) +{ + // + // The order in which the parameters are set is important + // since setCurvature might have to adjust the others. + // + mH = (h>=0) ? 1 : -1; // Default is: positive particle + // positive field + mOrigin = o; + setDipAngle(dip); + setPhase(phase); + + // + // Check for singularity and correct for negative curvature. + // May change mH and mPhase. Must therefore be set last. + // + setCurvature(c); + + // + // For the case B=0, h is ill defined. In the following we + // always assume h = +1. Since phase = psi - h * pi/2 + // we have to correct the phase in case h = -1. + // This assumes that the user uses the same h for phase + // as the one he passed to the constructor. + // + if (mSingularity && mH == -1) { + mH = +1; + setPhase(mPhase-M_PI); + } +} + +void Helix::setCurvature(double val) +{ + if (val < 0) { + mCurvature = -val; + mH = -mH; + setPhase(mPhase+M_PI); + } + else + mCurvature = val; + + if (fabs(mCurvature) <= std::numeric_limits::epsilon()) + mSingularity = true; // straight line + else + mSingularity = false; // curved +} + +void Helix::setPhase(double val) +{ + mPhase = val; + mCosPhase = cos(mPhase); + mSinPhase = sin(mPhase); + if (fabs(mPhase) > M_PI) + mPhase = atan2(mSinPhase, mCosPhase); // force range [-pi,pi] +} + +void Helix::setDipAngle(double val) +{ + mDipAngle = val; + mCosDipAngle = cos(mDipAngle); + mSinDipAngle = sin(mDipAngle); +} + +double Helix::xcenter() const +{ + if (mSingularity) + return 0; + else + return mOrigin.x-mCosPhase/mCurvature; +} + +double Helix::ycenter() const +{ + if (mSingularity) + return 0; + else + return mOrigin.y-mSinPhase/mCurvature; +} + +double Helix::fudgePathLength(const edm4hep::Vector3f& p) const +{ + double s; + double dx = p.x-mOrigin.x; + double dy = p.y-mOrigin.y; + + if (mSingularity) { + s = (dy*mCosPhase - dx*mSinPhase)/mCosDipAngle; + } + else { + s = atan2(dy*mCosPhase - dx*mSinPhase, + 1/mCurvature + dx*mCosPhase+dy*mSinPhase)/ + (mH*mCurvature*mCosDipAngle); + } + return s; +} + +double Helix::distance(const edm4hep::Vector3f& p, bool scanPeriods) const +{ + return edm4hep::utils::magnitude(this->at(pathLength(p,scanPeriods))-p); +} + +double Helix::pathLength(const edm4hep::Vector3f& p, bool scanPeriods) const +{ + // + // Returns the path length at the distance of closest + // approach between the helix and point p. + // For the case of B=0 (straight line) the path length + // can be calculated analytically. For B>0 there is + // unfortunately no easy solution to the problem. + // Here we use the Newton method to find the root of the + // referring equation. The 'fudgePathLength' serves + // as a starting value. + // + double s; + double dx = p.x-mOrigin.x; + double dy = p.y-mOrigin.y; + double dz = p.z-mOrigin.z; + + if (mSingularity) { + s = mCosDipAngle*(mCosPhase*dy-mSinPhase*dx) + + mSinDipAngle*dz; + } + else { // + const double MaxPrecisionNeeded = edm4eic::unit::um; + const int MaxIterations = 100; + + // + // The math is taken from Maple with C(expr,optimized) and + // some hand-editing. It is not very nice but efficient. + // + double t34 = mCurvature*mCosDipAngle*mCosDipAngle; + double t41 = mSinDipAngle*mSinDipAngle; + double t6, t7, t11, t12, t19; + + // + // Get a first guess by using the dca in 2D. Since + // in some extreme cases we might be off by n periods + // we add (subtract) periods in case we get any closer. + // + s = fudgePathLength(p); + + if (scanPeriods) { + double ds = period(); + int j, jmin = 0; + double d, dmin = edm4hep::utils::magnitude(at(s) - p); + for(j=1; j::max(); + else + return fabs(2*M_PI/(mH*mCurvature*mCosDipAngle)); +} + +std::pair Helix::pathLength(double r) const +{ + std::pair value; + std::pair VALUE(999999999.,999999999.); + // + // The math is taken from Maple with C(expr,optimized) and + // some hand-editing. It is not very nice but efficient. + // 'first' is the smallest of the two solutions (may be negative) + // 'second' is the other. + // + if (mSingularity) { + double t1 = mCosDipAngle*(mOrigin.x*mSinPhase-mOrigin.y*mCosPhase); + double t12 = mOrigin.y*mOrigin.y; + double t13 = mCosPhase*mCosPhase; + double t15 = r*r; + double t16 = mOrigin.x*mOrigin.x; + double t20 = -mCosDipAngle*mCosDipAngle*(2.0*mOrigin.x*mSinPhase*mOrigin.y*mCosPhase + + t12-t12*t13-t15+t13*t16); + if (t20<0.) return VALUE; + t20 = ::sqrt(t20); + value.first = (t1-t20)/(mCosDipAngle*mCosDipAngle); + value.second = (t1+t20)/(mCosDipAngle*mCosDipAngle); + } + else { + double t1 = mOrigin.y*mCurvature; + double t2 = mSinPhase; + double t3 = mCurvature*mCurvature; + double t4 = mOrigin.y*t2; + double t5 = mCosPhase; + double t6 = mOrigin.x*t5; + double t8 = mOrigin.x*mOrigin.x; + double t11 = mOrigin.y*mOrigin.y; + double t14 = r*r; + double t15 = t14*mCurvature; + double t17 = t8*t8; + double t19 = t11*t11; + double t21 = t11*t3; + double t23 = t5*t5; + double t32 = t14*t14; + double t35 = t14*t3; + double t38 = 8.0*t4*t6 - 4.0*t1*t2*t8 - 4.0*t11*mCurvature*t6 + + 4.0*t15*t6 + t17*t3 + t19*t3 + 2.0*t21*t8 + 4.0*t8*t23 - + 4.0*t8*mOrigin.x*mCurvature*t5 - 4.0*t11*t23 - + 4.0*t11*mOrigin.y*mCurvature*t2 + 4.0*t11 - 4.0*t14 + + t32*t3 + 4.0*t15*t4 - 2.0*t35*t11 - 2.0*t35*t8; + double t40 = (-t3*t38); + if (t40<0.) return VALUE; + t40 = ::sqrt(t40); + + double t43 = mOrigin.x*mCurvature; + double t45 = 2.0*t5 - t35 + t21 + 2.0 - 2.0*t1*t2 -2.0*t43 - 2.0*t43*t5 + t8*t3; + double t46 = mH*mCosDipAngle*mCurvature; + + value.first = (-mPhase + 2.0*atan((-2.0*t1 + 2.0*t2 + t40)/t45))/t46; + value.second = -(mPhase + 2.0*atan((2.0*t1 - 2.0*t2 + t40)/t45))/t46; + + // + // Solution can be off by +/- one period, select smallest + // + double p = period(); + if (! std::isnan(value.first)) { + if (fabs(value.first-p) < fabs(value.first)) value.first = value.first-p; + else if (fabs(value.first+p) < fabs(value.first)) value.first = value.first+p; + } + if (! std::isnan(value.second)) { + if (fabs(value.second-p) < fabs(value.second)) value.second = value.second-p; + else if (fabs(value.second+p) < fabs(value.second)) value.second = value.second+p; + } + } + if (value.first > value.second) + std::swap(value.first,value.second); + return(value); +} + +std::pair Helix::pathLength(double r, double x, double y) +{ + double x0 = mOrigin.x; + double y0 = mOrigin.y; + mOrigin.x = x0-x; + mOrigin.y = y0-y; + std::pair result = this->pathLength(r); + mOrigin.x = x0; + mOrigin.y = y0; + return result; +} + +double Helix::pathLength(const edm4hep::Vector3f& r, + const edm4hep::Vector3f& n) const +{ + // + // Vector 'r' defines the position of the center and + // vector 'n' the normal vector of the plane. + // For a straight line there is a simple analytical + // solution. For curvatures > 0 the root is determined + // by Newton method. In case no valid s can be found + // the max. largest value for s is returned. + // + double s; + + if (mSingularity) { + double t = n.z*mSinDipAngle + + n.y*mCosDipAngle*mCosPhase - + n.x*mCosDipAngle*mSinPhase; + if (t == 0) + s = NoSolution; + else + s = ((r - mOrigin)*n)/t; + } + else { + const double MaxPrecisionNeeded = edm4eic::unit::um; + const int MaxIterations = 20; + + double A = mCurvature*((mOrigin - r)*n) - + n.x*mCosPhase - + n.y*mSinPhase; + double t = mH*mCurvature*mCosDipAngle; + double u = n.z*mCurvature*mSinDipAngle; + + double a, f, fp; + double sOld = s = 0; + double shiftOld = 0; + double shift; +// (cos(angMax)-1)/angMax = 0.1 + const double angMax = 0.21; + double deltas = fabs(angMax/(mCurvature*mCosDipAngle)); +// dampingFactor = exp(-0.5); +// double dampingFactor = 0.60653; + int i; + + for (i=0; i +Helix::pathLengths(const Helix& h, double minStepSize, double minRange) const +{ + // + // Cannot handle case where one is a helix + // and the other one is a straight line. + // + if (mSingularity != h.mSingularity) + return std::pair(NoSolution, NoSolution); + + double s1, s2; + + if (mSingularity) { + // + // Analytic solution + // + edm4hep::Vector3f dv = h.mOrigin - mOrigin; + edm4hep::Vector3f a(-mCosDipAngle*mSinPhase, + mCosDipAngle*mCosPhase, + mSinDipAngle); + edm4hep::Vector3f b(-h.mCosDipAngle*h.mSinPhase, + h.mCosDipAngle*h.mCosPhase, + h.mSinDipAngle); + double ab = a*b; + double g = dv*a; + double k = dv*b; + s2 = (k-ab*g)/(ab*ab-1.); + s1 = g+s2*ab; + return std::pair(s1, s2); + } + else { + // + // First step: get dca in the xy-plane as start value + // + double dx = h.xcenter() - xcenter(); + double dy = h.ycenter() - ycenter(); + double dd = ::sqrt(dx*dx + dy*dy); + double r1 = 1/curvature(); + double r2 = 1/h.curvature(); + + double cosAlpha = (r1*r1 + dd*dd - r2*r2)/(2*r1*dd); + + double s; + double x, y; + if (fabs(cosAlpha) < 1) { // two solutions + double sinAlpha = sin(acos(cosAlpha)); + x = xcenter() + r1*(cosAlpha*dx - sinAlpha*dy)/dd; + y = ycenter() + r1*(sinAlpha*dx + cosAlpha*dy)/dd; + s = pathLength(x, y); + x = xcenter() + r1*(cosAlpha*dx + sinAlpha*dy)/dd; + y = ycenter() + r1*(cosAlpha*dy - sinAlpha*dx)/dd; + double a = pathLength(x, y); + if (h.distance(at(a)) < h.distance(at(s))) s = a; + } + else { // no intersection (or exactly one) + int rsign = ((r2-r1) > dd ? -1 : 1); // set -1 when *this* helix is + // completely contained in the other + x = xcenter() + rsign*r1*dx/dd; + y = ycenter() + rsign*r1*dy/dd; + s = pathLength(x, y); + } + + // + // Second step: scan in decreasing intervals around seed 's' + // minRange and minStepSize are passed as arguments to the method. + // They have default values defined in the header file. + // + double dmin = h.distance(at(s)); + double range = std::max(2*dmin, minRange); + double ds = range/10; + double slast=-999999, ss, d; + s1 = s - range/2.; + s2 = s + range/2.; + + while (ds > minStepSize) { + for (ss=s1; ss(s, h.pathLength(at(s))); + } +} + + +void Helix::moveOrigin(double s) +{ + if (mSingularity) + mOrigin = at(s); + else { + edm4hep::Vector3f newOrigin = at(s); + double newPhase = atan2(newOrigin.y - ycenter(), + newOrigin.x - xcenter()); + mOrigin = newOrigin; + setPhase(newPhase); + } +} +/* +int operator== (const Helix& a, const Helix& b) +{ + // + // Checks for numerical identity only ! + // + return (a.origin() == b.origin() && + a.dipAngle() == b.dipAngle() && + a.curvature() == b.curvature() && + a.phase() == b.phase() && + a.h() == b.h()); +} + +int operator!= (const Helix& a, const Helix& b) {return !(a == b);} +*/ +//std::ostream& operator<<(std::ostream& os, const Helix& h) +void Helix::Print() const +{ + std::cout << "(" + << "curvature = " << mCurvature << ", " + << "dip angle = " << mDipAngle << ", " + << "phase = " << mPhase << ", " + << "h = " << mH << ", " + << "origin = " << mOrigin.x << " " << mOrigin.y << " " << mOrigin.z << ")" + << std::endl; +} + +edm4hep::Vector3f Helix::momentum(double B) const +{ + if (mSingularity) + return(edm4hep::Vector3f(0,0,0)); + else { + double pt = edm4eic::unit::GeV*fabs(dd4hep::c_light*dd4hep::nanosecond/dd4hep::meter*B/dd4hep::tesla)/(fabs(mCurvature)*dd4hep::meter); + + return (edm4hep::Vector3f(pt*cos(mPhase+mH*M_PI/2), // pos part pos field + pt*sin(mPhase+mH*M_PI/2), + pt*tan(mDipAngle))); + } +} + +edm4hep::Vector3f Helix::momentumAt(double S, double B) const +{ + // Obtain phase-shifted momentum from phase-shift of origin + Helix tmp(*this); + tmp.moveOrigin(S); + return tmp.momentum(B); +} + +int Helix::charge(double B) const +{ + return (B > 0 ? -mH : mH); +} + +double Helix::geometricSignedDistance(double x, double y) +{ + // Geometric signed distance + double thePath = this->pathLength(x,y); + edm4hep::Vector3f DCA2dPosition = this->at(thePath); + DCA2dPosition.z = 0; + edm4hep::Vector3f position(x,y,0); + edm4hep::Vector3f DCAVec = (DCA2dPosition-position); + edm4hep::Vector3f momVec; + // Deal with straight tracks + if (this->mSingularity) { + momVec = this->at(1)- this->at(0); + momVec.z = 0; + } + else { + momVec = this->momentumAt(thePath,1./dd4hep::tesla); // Don't care about Bmag. Helicity is what matters. + momVec.z = 0; + } + + double cross = DCAVec.x*momVec.y - DCAVec.y*momVec.x; + double theSign = (cross>=0) ? 1. : -1.; + return theSign*edm4hep::utils::magnitudeTransverse(DCAVec); +} + +double Helix::curvatureSignedDistance(double x, double y) +{ + // Protect against mH = 0 or zero field + if (this->mSingularity || abs(this->mH)<=0) { + return (this->geometricSignedDistance(x,y)); + } + else { + return (this->geometricSignedDistance(x,y))/(this->mH); + } + +} + +double Helix::geometricSignedDistance(const edm4hep::Vector3f& pos) +{ + double sdca2d = this->geometricSignedDistance(pos.x,pos.y); + double theSign = (sdca2d>=0) ? 1. : -1.; + return (this->distance(pos))*theSign; +} + +double Helix::curvatureSignedDistance(const edm4hep::Vector3f& pos) +{ + double sdca2d = this->curvatureSignedDistance(pos.x,pos.y); + double theSign = (sdca2d>=0) ? 1. : -1.; + return (this->distance(pos))*theSign; +} + +} // namespace eicrecon \ No newline at end of file diff --git a/src/algorithms/reco/Helix.h b/src/algorithms/reco/Helix.h new file mode 100644 index 0000000000..37a281d6b5 --- /dev/null +++ b/src/algorithms/reco/Helix.h @@ -0,0 +1,238 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Xin Dong, Rongrong Ma + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace eicrecon { + +class Helix { + bool mSingularity; // true for straight line case (B=0) + edm4hep::Vector3f mOrigin; + double mDipAngle; + double mCurvature; + double mPhase; + int mH; // -sign(q*B); + + double mCosDipAngle; + double mSinDipAngle; + double mCosPhase; + double mSinPhase; + +public: + /// curvature, dip angle, phase, origin, h + Helix(const double c, const double dip, const double phase, const edm4hep::Vector3f& o, const int h=-1); + + /// momentum, origin, b_field, charge + Helix(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q); + + /// ReconstructParticle, b field + Helix(const edm4eic::ReconstructedParticle& p, const double b_field); + + ~Helix() = default; + + double dipAngle() const; + double curvature() const; /// 1/R in xy-plane + double phase() const; /// aziumth in xy-plane measured from ring center + double xcenter() const; /// x-center of circle in xy-plane + double ycenter() const; /// y-center of circle in xy-plane + int h() const; /// -sign(q*B); + + const edm4hep::Vector3f& origin() const; /// starting point + + void setParameters(double c, double dip, double phase, const edm4hep::Vector3f& o, int h); + + void setParameters(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q); + + /// edm4eic::TrackParameters, b field + void setParameters(const edm4eic::TrackParameters& trk, const double b_field); + + /// coordinates of helix at point s + double x(double s) const; + double y(double s) const; + double z(double s) const; + + edm4hep::Vector3f at(double s) const; + + /// pointing vector of helix at point s + double cx(double s) const; + double cy(double s) const; + double cz(double s = 0) const; + + edm4hep::Vector3f cat(double s) const; + + /// returns period length of helix + double period() const; + + /// path length at given r (cylindrical r) + std::pair pathLength(double r) const; + + /// path length at given r (cylindrical r, cylinder axis at x,y) + std::pair pathLength(double r, double x, double y); + + /// path length at distance of closest approach to a given point + double pathLength(const edm4hep::Vector3f& p, bool scanPeriods = true) const; + + /// path length at intersection with plane + double pathLength(const edm4hep::Vector3f& r, + const edm4hep::Vector3f& n) const; + + /// path length at distance of closest approach in the xy-plane to a given point + double pathLength(double x, double y) const; + + /// path lengths at dca between two helices + std::pair pathLengths(const Helix&, + double minStepSize = 10*edm4eic::unit::um, + double minRange = 10*edm4eic::unit::cm) const; + + /// minimal distance between point and helix + double distance(const edm4hep::Vector3f& p, bool scanPeriods = true) const; + + /// checks for valid parametrization + bool valid(double world = 1.e+5) const {return !bad(world);} + int bad(double world = 1.e+5) const; + + /// move the origin along the helix to s which becomes then s=0 + virtual void moveOrigin(double s); + + static const double NoSolution; + + + void setCurvature(double); /// performs also various checks + void setPhase(double); + void setDipAngle(double); + + /// value of S where distance in x-y plane is minimal + double fudgePathLength(const edm4hep::Vector3f&) const; + + // Requires: signed Magnetic Field + edm4hep::Vector3f momentum(double) const; // returns the momentum at origin + edm4hep::Vector3f momentumAt(double, double) const; // returns momemtum at S + int charge(double) const; // returns charge of particle + // 2d DCA to x,y point signed relative to curvature + double curvatureSignedDistance(double x, double y) ; + // 2d DCA to x,y point signed relative to rotation + double geometricSignedDistance(double x, double y) ; + // 3d DCA to 3d point signed relative to curvature + double curvatureSignedDistance(const edm4hep::Vector3f&) ; + // 3d DCA to 3d point signed relative to rotation + double geometricSignedDistance(const edm4hep::Vector3f&) ; + + // + void Print() const; + +}; // end class Helix + +// +// Non-member functions +// +//int operator== (const Helix&, const Helix&); +//int operator!= (const Helix&, const Helix&); +//std::ostream& operator<<(std::ostream&, const Helix&); + +// +// Inline functions +// +inline int Helix::h() const {return mH;} + +inline double Helix::dipAngle() const {return mDipAngle;} + +inline double Helix::curvature() const {return mCurvature;} + +inline double Helix::phase() const {return mPhase;} + +inline double Helix::x(double s) const +{ + if (mSingularity) + return mOrigin.x - s*mCosDipAngle*mSinPhase; + else + return mOrigin.x + (cos(mPhase + s*mH*mCurvature*mCosDipAngle)-mCosPhase)/mCurvature; +} + +inline double Helix::y(double s) const +{ + if (mSingularity) + return mOrigin.y + s*mCosDipAngle*mCosPhase; + else + return mOrigin.y + (sin(mPhase + s*mH*mCurvature*mCosDipAngle)-mSinPhase)/mCurvature; +} + +inline double Helix::z(double s) const +{ + return mOrigin.z + s*mSinDipAngle; +} + +inline double Helix::cx(double s) const +{ + if (mSingularity) + return -mCosDipAngle*mSinPhase; + else + return -sin(mPhase + s*mH*mCurvature*mCosDipAngle)*mH*mCosDipAngle; +} + +inline double Helix::cy(double s) const +{ + if (mSingularity) + return mCosDipAngle*mCosPhase; + else + return cos(mPhase + s*mH*mCurvature*mCosDipAngle)*mH*mCosDipAngle; +} + +inline double Helix::cz(double /* s */) const +{ + return mSinDipAngle; +} + +inline const edm4hep::Vector3f& Helix::origin() const {return mOrigin;} + +inline edm4hep::Vector3f Helix::at(double s) const +{ + return edm4hep::Vector3f(x(s), y(s), z(s)); +} + +inline edm4hep::Vector3f Helix::cat(double s) const +{ + return edm4hep::Vector3f(cx(s), cy(s), cz(s)); +} + +inline double Helix::pathLength(double X, double Y) const +{ + return fudgePathLength(edm4hep::Vector3f(X, Y, 0)); +} +inline int Helix::bad(double WorldSize) const +{ + +// int ierr; + if (!::finite(mDipAngle )) return 11; + if (!::finite(mCurvature )) return 12; + +// ierr = mOrigin.bad(WorldSize); +// if (ierr) return 3+ierr*100; + + if (::fabs(mDipAngle) >1.58) return 21; + double qwe = ::fabs(::fabs(mDipAngle)-M_PI/2); + if (qwe < 1./WorldSize ) return 31; + + if (::fabs(mCurvature) > WorldSize) return 22; + if (mCurvature < 0 ) return 32; + + if (abs(mH) != 1 ) return 24; + + return 0; +} + +} diff --git a/src/algorithms/reco/SecondaryVertices.cc b/src/algorithms/reco/SecondaryVertices.cc new file mode 100644 index 0000000000..50b377aa05 --- /dev/null +++ b/src/algorithms/reco/SecondaryVertices.cc @@ -0,0 +1,125 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Daniel Brandenburg, Xin Dong + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "algorithms/reco/SecondaryVertices.h" +#include "algorithms/reco/Helix.h" +#include "algorithms/reco/SecondaryVerticesConfig.h" +#include "services/particle/ParticleSvc.h" + +namespace eicrecon { + +/** + * @brief Initialize the SecondaryVertices Algorithm + * + */ +void SecondaryVertices::init() {} + +/** + * @brief Produce a list of secondary vertex candidates + * + * @param rcvtx - input collection of all vertex candidates + * @return edm4eic::VertexCollection + */ +void SecondaryVertices::process(const SecondaryVertices::Input& input, + const SecondaryVertices::Output& output) const { + const auto [rcvtx, rcparts] = input; + auto [out_secondary_vertices] = output; + + auto& particleSvc = algorithms::ParticleSvc::instance(); +// edm4hep::Vector3f pVtxPos; +// for(const auto& v : primVtx) + const auto pVtxPos4f = (*rcvtx)[0].getPosition(); + // convert to cm + edm4hep::Vector3f pVtxPos(pVtxPos4f.x*edm4eic::unit::mm/edm4eic::unit::cm, + pVtxPos4f.y*edm4eic::unit::mm/edm4eic::unit::cm, + pVtxPos4f.z*edm4eic::unit::mm/edm4eic::unit::cm); + info("\t Primary vertex = ({},{},{})cm \t b field = {} tesla", pVtxPos.x, pVtxPos.y, pVtxPos.z, m_cfg.b_field/dd4hep::tesla); + + std::vector pi_index; + std::vector k_index; + std::vector p_index; + for (unsigned int i = 0; const auto& p : *rcparts) { + const auto pdg = p.getPDG(); + if(abs(pdg) == 211) pi_index.push_back(i); + if(abs(pdg) == 321) k_index.push_back(i); + if(abs(pdg) == 2212) p_index.push_back(i); + ++i; + } + + info("\t Array sizes: pions = {}, kaons = {}, protons = {}", pi_index.size(), k_index.size(), p_index.size()); + + for (unsigned int i1 = 0; i1 < pi_index.size(); ++i1) { + for (unsigned int i2 = i1 + 1; i2 < pi_index.size(); ++i2) { + const auto& p1 = (*rcparts)[i1]; + const auto& p2 = (*rcparts)[i2]; + + if (p1.getCharge() + p2.getCharge() != 0) continue; + + Helix h1obj(p1, m_cfg.b_field); Helix& h1 = h1obj; + Helix h2obj(p2, m_cfg.b_field); Helix& h2 = h2obj; + + // Helix function uses cm unit + double dca1 = h1.distance(pVtxPos) * edm4eic::unit::cm; + double dca2 = h2.distance(pVtxPos) * edm4eic::unit::cm; + debug("\t dca1 = {}, dca2 = {}", dca1, dca2); + if( dca1 < m_cfg.minDca1 || dca2 < m_cfg.minDca2 ) continue; + + std::pair const ss = h1.pathLengths(h2); + edm4hep::Vector3f h1AtDcaTo2 = h1.at(ss.first); + edm4hep::Vector3f h2AtDcaTo1 = h2.at(ss.second); + + double dca12 = edm4hep::utils::magnitude(h1AtDcaTo2 - h2AtDcaTo1) * edm4eic::unit::cm; + if( dca12 > m_cfg.maxDca12 ) continue; + edm4hep::Vector3f pairPos = 0.5*(h1AtDcaTo2 + h2AtDcaTo1); + + edm4hep::Vector3f h1MomAtDca = h1.momentumAt(ss.first, m_cfg.b_field); + edm4hep::Vector3f h2MomAtDca = h2.momentumAt(ss.second, m_cfg.b_field); + edm4hep::Vector3f pairMom = h1MomAtDca + h2MomAtDca; + + float e1 = std::hypot(edm4hep::utils::magnitude(h1MomAtDca), particleSvc.particle(211).mass); + float e2 = std::hypot(edm4hep::utils::magnitude(h2MomAtDca), particleSvc.particle(211).mass); + float pairE = e1+e2; + + edm4hep::Vector4f h1FourMom(h1MomAtDca.x, h1MomAtDca.y, h1MomAtDca.z, e1); + edm4hep::Vector4f h2FourMom(h2MomAtDca.x, h2MomAtDca.y, h2MomAtDca.z, e2); + + double m_inv = std::hypot(pairE, -edm4hep::utils::magnitude(pairMom)); + double angle = edm4hep::utils::angleBetween(pairMom, pairPos - pVtxPos); + if(cos(angle) < m_cfg.minCostheta ) continue; + + double beta = edm4hep::utils::magnitude(pairMom)/pairE; + double time = edm4hep::utils::magnitude(pairPos - pVtxPos)/(beta*dd4hep::c_light); + auto v0 = out_secondary_vertices->create(); + v0.setType(2); // 2 for secondary + v0.setPosition({(float)(pairPos.x * edm4eic::unit::cm / edm4eic::unit::mm), + (float)(pairPos.y * edm4eic::unit::cm / edm4eic::unit::mm), + (float)(pairPos.z * edm4eic::unit::cm / edm4eic::unit::mm), + (float)time}); + v0.addToAssociatedParticles(p1); + v0.addToAssociatedParticles(p2); + + info("One secondary vertex found at (x,y,z) = ({}, {}, {}) mm.", + pairPos.x * edm4eic::unit::cm / edm4eic::unit::mm, + pairPos.y * edm4eic::unit::cm / edm4eic::unit::mm, + pairPos.x * edm4eic::unit::cm / edm4eic::unit::mm); + + } // end i2 + } // end i1 + +} // end process + +} // namespace eicrecon diff --git a/src/algorithms/reco/SecondaryVertices.h b/src/algorithms/reco/SecondaryVertices.h new file mode 100644 index 0000000000..35bb05420f --- /dev/null +++ b/src/algorithms/reco/SecondaryVertices.h @@ -0,0 +1,37 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Daniel Brandenburg, Xin Dong + +#pragma once + +#include +#include +#include +#include // for basic_string +#include // for string_view + +#include "algorithms/interfaces/WithPodConfig.h" +#include "algorithms/reco/SecondaryVerticesConfig.h" + +namespace eicrecon { + +using SecondaryVerticesAlgorithm = algorithms::Algorithm< + algorithms::Input, + algorithms::Output>; + +class SecondaryVertices : public SecondaryVerticesAlgorithm, + public WithPodConfig { + +public: + SecondaryVertices(std::string_view name) + : SecondaryVerticesAlgorithm{name, + {"inputVertices", "inputParticles"}, + {"outputSecondaryVertices"}, + "Reconstruct secondary vertices in SecondaryVertices collection"} {} + + void init() final; + void process(const Input&, const Output&) const final; + +private: + SecondaryVerticesConfig m_cfg; +}; +} // namespace eicrecon diff --git a/src/algorithms/reco/SecondaryVerticesConfig.h b/src/algorithms/reco/SecondaryVerticesConfig.h new file mode 100644 index 0000000000..daedee6fc6 --- /dev/null +++ b/src/algorithms/reco/SecondaryVerticesConfig.h @@ -0,0 +1,22 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Xin Dong + +#pragma once + +#include +#include +#include + +namespace eicrecon { + +struct SecondaryVerticesConfig { + + float b_field = -1.7*dd4hep::tesla; + float minDca1 = 0.03*edm4eic::unit::mm; // mm, daughter1 to pVtx + float minDca2 = 0.03*edm4eic::unit::mm; // mm, daughter2 to pVtx + float maxDca12 = 0.3*edm4eic::unit::mm; // mm, dca between daughter 1 and 2 + float maxDca = 0.3*edm4eic::unit::mm; // mm, dca of V0 to pVtx + float minCostheta = 0.95; // costheta, theta: angle of V0 decay direction and momentum +}; + +} // namespace eicrecon diff --git a/src/factories/reco/SecondaryVertices_factory.h b/src/factories/reco/SecondaryVertices_factory.h new file mode 100644 index 0000000000..6f080f0832 --- /dev/null +++ b/src/factories/reco/SecondaryVertices_factory.h @@ -0,0 +1,50 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Xin Dong + +#pragma once + +#include +#include +#include +#include +#include + +#include "algorithms/reco/SecondaryVertices.h" +#include "algorithms/reco/SecondaryVerticesConfig.h" +#include "services/algorithms_init/AlgorithmsInit_service.h" +#include "extensions/jana/JOmniFactory.h" + +namespace eicrecon { + +class SecondaryVertices_factory + : public JOmniFactory { + +public: + using AlgoT = eicrecon::SecondaryVertices; + +private: + std::unique_ptr m_algo; + + PodioInput m_rc_vertices_input{this}; + PodioInput m_rc_parts_input{this}; + + // Declare outputs + PodioOutput m_secondary_vertices_output{this}; + + Service m_algorithmsInit{this}; + +public: + void Configure() { + m_algo = std::make_unique(GetPrefix()); + m_algo->level((algorithms::LogLevel)logger()->level()); + + m_algo->applyConfig(config()); + m_algo->init(); + } + + void Process(int32_t /* run_number */, uint64_t /* event_number */) { + m_algo->process({m_rc_vertices_input(), m_rc_parts_input()}, {m_secondary_vertices_output().get()}); + } +}; + +} // namespace eicrecon diff --git a/src/global/reco/reco.cc b/src/global/reco/reco.cc index f2af731a7d..9b7b0b2c60 100644 --- a/src/global/reco/reco.cc +++ b/src/global/reco/reco.cc @@ -38,6 +38,7 @@ #include "factories/reco/MC2ReconstructedParticle_factory.h" #include "factories/reco/MatchClusters_factory.h" #include "factories/reco/PrimaryVertices_factory.h" +#include "factories/reco/SecondaryVertices_factory.h" #include "factories/reco/ReconstructedElectrons_factory.h" #include "factories/reco/ScatteredElectronsEMinusPz_factory.h" #include "factories/reco/ScatteredElectronsTruth_factory.h" @@ -77,7 +78,6 @@ void InitPlugin(JApplication* app) { "EcalEndcapPClusterAssociations"}, {"EcalClusterAssociations"}, app)); - // Create ReconstructedParticles app->Add(new JOmniFactoryGeneratorT( "ReconstructedParticlesWithAssoc", { @@ -274,5 +274,9 @@ void InitPlugin(JApplication* app) { app->Add(new JOmniFactoryGeneratorT( "PrimaryVertices", {"CentralTrackVertices"}, {"PrimaryVertices"}, {}, app)); + + app->Add(new JOmniFactoryGeneratorT( + "SecondaryVertices", {"PrimaryVertices","ReconstructedParticles"}, {"SecondaryVertices"}, {}, app)); + } } // extern "C" diff --git a/src/services/io/podio/JEventProcessorPODIO.cc b/src/services/io/podio/JEventProcessorPODIO.cc index e77a77f3f5..32daecd5bf 100644 --- a/src/services/io/podio/JEventProcessorPODIO.cc +++ b/src/services/io/podio/JEventProcessorPODIO.cc @@ -258,6 +258,7 @@ JEventProcessorPODIO::JEventProcessorPODIO() { "ScatteredElectronsTruth", "ScatteredElectronsEMinusPz", "PrimaryVertices", + "SecondaryVertices", "BarrelClusters", "HadronicFinalState", From f04e71a2863c8f585dc77a81ac7c758ad8f7caf9 Mon Sep 17 00:00:00 2001 From: Xin Dong Date: Thu, 9 Oct 2025 14:03:04 -0400 Subject: [PATCH 02/17] rename SecondaryVertices to SecondaryVerticesHelix --- ...econdaryVertices.cc => SecondaryVerticesHelix.cc} | 12 ++++++------ ...{SecondaryVertices.h => SecondaryVerticesHelix.h} | 12 ++++++------ ...rticesConfig.h => SecondaryVerticesHelixConfig.h} | 2 +- ...es_factory.h => SecondaryVerticesHelix_factory.h} | 10 +++++----- src/global/reco/reco.cc | 6 +++--- src/services/io/podio/JEventProcessorPODIO.cc | 2 +- 6 files changed, 22 insertions(+), 22 deletions(-) rename src/algorithms/reco/{SecondaryVertices.cc => SecondaryVerticesHelix.cc} (92%) rename src/algorithms/reco/{SecondaryVertices.h => SecondaryVerticesHelix.h} (73%) rename src/algorithms/reco/{SecondaryVerticesConfig.h => SecondaryVerticesHelixConfig.h} (94%) rename src/factories/reco/{SecondaryVertices_factory.h => SecondaryVerticesHelix_factory.h} (78%) diff --git a/src/algorithms/reco/SecondaryVertices.cc b/src/algorithms/reco/SecondaryVerticesHelix.cc similarity index 92% rename from src/algorithms/reco/SecondaryVertices.cc rename to src/algorithms/reco/SecondaryVerticesHelix.cc index 50b377aa05..0bc8d1a77c 100644 --- a/src/algorithms/reco/SecondaryVertices.cc +++ b/src/algorithms/reco/SecondaryVerticesHelix.cc @@ -15,18 +15,18 @@ #include #include -#include "algorithms/reco/SecondaryVertices.h" +#include "algorithms/reco/SecondaryVerticesHelix.h" #include "algorithms/reco/Helix.h" -#include "algorithms/reco/SecondaryVerticesConfig.h" +#include "algorithms/reco/SecondaryVerticesHelixConfig.h" #include "services/particle/ParticleSvc.h" namespace eicrecon { /** - * @brief Initialize the SecondaryVertices Algorithm + * @brief Initialize the SecondaryVerticesHelix Algorithm * */ -void SecondaryVertices::init() {} +void SecondaryVerticesHelix::init() {} /** * @brief Produce a list of secondary vertex candidates @@ -34,8 +34,8 @@ void SecondaryVertices::init() {} * @param rcvtx - input collection of all vertex candidates * @return edm4eic::VertexCollection */ -void SecondaryVertices::process(const SecondaryVertices::Input& input, - const SecondaryVertices::Output& output) const { +void SecondaryVerticesHelix::process(const SecondaryVerticesHelix::Input& input, + const SecondaryVerticesHelix::Output& output) const { const auto [rcvtx, rcparts] = input; auto [out_secondary_vertices] = output; diff --git a/src/algorithms/reco/SecondaryVertices.h b/src/algorithms/reco/SecondaryVerticesHelix.h similarity index 73% rename from src/algorithms/reco/SecondaryVertices.h rename to src/algorithms/reco/SecondaryVerticesHelix.h index 35bb05420f..80e1622390 100644 --- a/src/algorithms/reco/SecondaryVertices.h +++ b/src/algorithms/reco/SecondaryVerticesHelix.h @@ -14,16 +14,16 @@ namespace eicrecon { -using SecondaryVerticesAlgorithm = algorithms::Algorithm< +using SecondaryVerticesHelixAlgorithm = algorithms::Algorithm< algorithms::Input, algorithms::Output>; -class SecondaryVertices : public SecondaryVerticesAlgorithm, - public WithPodConfig { +class SecondaryVerticesHelix : public SecondaryVerticesHelixAlgorithm, + public WithPodConfig { public: - SecondaryVertices(std::string_view name) - : SecondaryVerticesAlgorithm{name, + SecondaryVerticesHelix(std::string_view name) + : SecondaryVerticesHelixAlgorithm{name, {"inputVertices", "inputParticles"}, {"outputSecondaryVertices"}, "Reconstruct secondary vertices in SecondaryVertices collection"} {} @@ -32,6 +32,6 @@ class SecondaryVertices : public SecondaryVerticesAlgorithm, void process(const Input&, const Output&) const final; private: - SecondaryVerticesConfig m_cfg; + SecondaryVerticesHelixConfig m_cfg; }; } // namespace eicrecon diff --git a/src/algorithms/reco/SecondaryVerticesConfig.h b/src/algorithms/reco/SecondaryVerticesHelixConfig.h similarity index 94% rename from src/algorithms/reco/SecondaryVerticesConfig.h rename to src/algorithms/reco/SecondaryVerticesHelixConfig.h index daedee6fc6..00e47c6ae1 100644 --- a/src/algorithms/reco/SecondaryVerticesConfig.h +++ b/src/algorithms/reco/SecondaryVerticesHelixConfig.h @@ -9,7 +9,7 @@ namespace eicrecon { -struct SecondaryVerticesConfig { +struct SecondaryVerticesHelixConfig { float b_field = -1.7*dd4hep::tesla; float minDca1 = 0.03*edm4eic::unit::mm; // mm, daughter1 to pVtx diff --git a/src/factories/reco/SecondaryVertices_factory.h b/src/factories/reco/SecondaryVerticesHelix_factory.h similarity index 78% rename from src/factories/reco/SecondaryVertices_factory.h rename to src/factories/reco/SecondaryVerticesHelix_factory.h index 6f080f0832..17fa2df7b5 100644 --- a/src/factories/reco/SecondaryVertices_factory.h +++ b/src/factories/reco/SecondaryVerticesHelix_factory.h @@ -9,18 +9,18 @@ #include #include -#include "algorithms/reco/SecondaryVertices.h" -#include "algorithms/reco/SecondaryVerticesConfig.h" +#include "algorithms/reco/SecondaryVerticesHelix.h" +#include "algorithms/reco/SecondaryVerticesHelixConfig.h" #include "services/algorithms_init/AlgorithmsInit_service.h" #include "extensions/jana/JOmniFactory.h" namespace eicrecon { -class SecondaryVertices_factory - : public JOmniFactory { +class SecondaryVerticesHelix_factory + : public JOmniFactory { public: - using AlgoT = eicrecon::SecondaryVertices; + using AlgoT = eicrecon::SecondaryVerticesHelix; private: std::unique_ptr m_algo; diff --git a/src/global/reco/reco.cc b/src/global/reco/reco.cc index 9b7b0b2c60..a934afec99 100644 --- a/src/global/reco/reco.cc +++ b/src/global/reco/reco.cc @@ -38,7 +38,7 @@ #include "factories/reco/MC2ReconstructedParticle_factory.h" #include "factories/reco/MatchClusters_factory.h" #include "factories/reco/PrimaryVertices_factory.h" -#include "factories/reco/SecondaryVertices_factory.h" +#include "factories/reco/SecondaryVerticesHelix_factory.h" #include "factories/reco/ReconstructedElectrons_factory.h" #include "factories/reco/ScatteredElectronsEMinusPz_factory.h" #include "factories/reco/ScatteredElectronsTruth_factory.h" @@ -275,8 +275,8 @@ void InitPlugin(JApplication* app) { app->Add(new JOmniFactoryGeneratorT( "PrimaryVertices", {"CentralTrackVertices"}, {"PrimaryVertices"}, {}, app)); - app->Add(new JOmniFactoryGeneratorT( - "SecondaryVertices", {"PrimaryVertices","ReconstructedParticles"}, {"SecondaryVertices"}, {}, app)); + app->Add(new JOmniFactoryGeneratorT( + "SecondaryVerticesHelix", {"PrimaryVertices","ReconstructedParticles"}, {"SecondaryVerticesHelix"}, {}, app)); } } // extern "C" diff --git a/src/services/io/podio/JEventProcessorPODIO.cc b/src/services/io/podio/JEventProcessorPODIO.cc index 32daecd5bf..b3d7ecc3ca 100644 --- a/src/services/io/podio/JEventProcessorPODIO.cc +++ b/src/services/io/podio/JEventProcessorPODIO.cc @@ -258,7 +258,7 @@ JEventProcessorPODIO::JEventProcessorPODIO() { "ScatteredElectronsTruth", "ScatteredElectronsEMinusPz", "PrimaryVertices", - "SecondaryVertices", + "SecondaryVerticesHelix", "BarrelClusters", "HadronicFinalState", From e14e205a650567b1e5de28b6032139108f8c0668 Mon Sep 17 00:00:00 2001 From: Xin Dong Date: Thu, 9 Oct 2025 16:22:34 -0400 Subject: [PATCH 03/17] typo fix when renaming --- src/algorithms/reco/SecondaryVerticesHelix.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/reco/SecondaryVerticesHelix.h b/src/algorithms/reco/SecondaryVerticesHelix.h index 80e1622390..378cb45558 100644 --- a/src/algorithms/reco/SecondaryVerticesHelix.h +++ b/src/algorithms/reco/SecondaryVerticesHelix.h @@ -10,7 +10,7 @@ #include // for string_view #include "algorithms/interfaces/WithPodConfig.h" -#include "algorithms/reco/SecondaryVerticesConfig.h" +#include "algorithms/reco/SecondaryVerticesHelixConfig.h" namespace eicrecon { From ebe0657b49d0b7003e5ce7c029c918180156c1f6 Mon Sep 17 00:00:00 2001 From: Xin Dong Date: Fri, 17 Oct 2025 13:59:18 -0400 Subject: [PATCH 04/17] fix invariant mass calculation --- src/algorithms/reco/SecondaryVerticesHelix.cc | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/algorithms/reco/SecondaryVerticesHelix.cc b/src/algorithms/reco/SecondaryVerticesHelix.cc index 0bc8d1a77c..3e8dcf2e19 100644 --- a/src/algorithms/reco/SecondaryVerticesHelix.cc +++ b/src/algorithms/reco/SecondaryVerticesHelix.cc @@ -90,14 +90,13 @@ void SecondaryVerticesHelix::process(const SecondaryVerticesHelix::Input& input, edm4hep::Vector3f h2MomAtDca = h2.momentumAt(ss.second, m_cfg.b_field); edm4hep::Vector3f pairMom = h1MomAtDca + h2MomAtDca; - float e1 = std::hypot(edm4hep::utils::magnitude(h1MomAtDca), particleSvc.particle(211).mass); - float e2 = std::hypot(edm4hep::utils::magnitude(h2MomAtDca), particleSvc.particle(211).mass); - float pairE = e1+e2; + double e1 = std::hypot(edm4hep::utils::magnitude(h1MomAtDca), particleSvc.particle(211).mass); + double e2 = std::hypot(edm4hep::utils::magnitude(h2MomAtDca), particleSvc.particle(211).mass); + double pairE = e1+e2; + double pairP = edm4hep::utils::magnitude(pairMom); - edm4hep::Vector4f h1FourMom(h1MomAtDca.x, h1MomAtDca.y, h1MomAtDca.z, e1); - edm4hep::Vector4f h2FourMom(h2MomAtDca.x, h2MomAtDca.y, h2MomAtDca.z, e2); - - double m_inv = std::hypot(pairE, -edm4hep::utils::magnitude(pairMom)); + double m_inv2 = pairE*pairE - pairP*pairP; + double m_inv = (m_inv2>0) sqrt(m_inv2) : 0.; double angle = edm4hep::utils::angleBetween(pairMom, pairPos - pVtxPos); if(cos(angle) < m_cfg.minCostheta ) continue; From 9b81bab59d13d7c6661b087641524aa4fc903f20 Mon Sep 17 00:00:00 2001 From: Xin Dong Date: Fri, 17 Oct 2025 14:12:54 -0400 Subject: [PATCH 05/17] fix a typo --- src/algorithms/reco/SecondaryVerticesHelix.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/reco/SecondaryVerticesHelix.cc b/src/algorithms/reco/SecondaryVerticesHelix.cc index 3e8dcf2e19..b383be2a61 100644 --- a/src/algorithms/reco/SecondaryVerticesHelix.cc +++ b/src/algorithms/reco/SecondaryVerticesHelix.cc @@ -96,7 +96,7 @@ void SecondaryVerticesHelix::process(const SecondaryVerticesHelix::Input& input, double pairP = edm4hep::utils::magnitude(pairMom); double m_inv2 = pairE*pairE - pairP*pairP; - double m_inv = (m_inv2>0) sqrt(m_inv2) : 0.; + double m_inv = (m_inv2>0) ? sqrt(m_inv2) : 0.; double angle = edm4hep::utils::angleBetween(pairMom, pairPos - pVtxPos); if(cos(angle) < m_cfg.minCostheta ) continue; From 54a2cb24a74f9dd85f10de221519fddc66c4cef3 Mon Sep 17 00:00:00 2001 From: Xin Dong Date: Mon, 20 Oct 2025 17:10:14 -0400 Subject: [PATCH 06/17] Generalized the reconstruction for any two-track combination --- src/algorithms/reco/SecondaryVerticesHelix.cc | 64 +++++++++++-------- src/algorithms/reco/SecondaryVerticesHelix.h | 2 +- .../reco/SecondaryVerticesHelixConfig.h | 10 +-- 3 files changed, 44 insertions(+), 32 deletions(-) diff --git a/src/algorithms/reco/SecondaryVerticesHelix.cc b/src/algorithms/reco/SecondaryVerticesHelix.cc index b383be2a61..e1fc0d2ec7 100644 --- a/src/algorithms/reco/SecondaryVerticesHelix.cc +++ b/src/algorithms/reco/SecondaryVerticesHelix.cc @@ -40,49 +40,56 @@ void SecondaryVerticesHelix::process(const SecondaryVerticesHelix::Input& input, auto [out_secondary_vertices] = output; auto& particleSvc = algorithms::ParticleSvc::instance(); -// edm4hep::Vector3f pVtxPos; -// for(const auto& v : primVtx) + + if((*rcvtx).size()==0) { + info(" No primary vertex in this event! Skip secondary vertex finder!"); + return; + } const auto pVtxPos4f = (*rcvtx)[0].getPosition(); // convert to cm edm4hep::Vector3f pVtxPos(pVtxPos4f.x*edm4eic::unit::mm/edm4eic::unit::cm, pVtxPos4f.y*edm4eic::unit::mm/edm4eic::unit::cm, pVtxPos4f.z*edm4eic::unit::mm/edm4eic::unit::cm); info("\t Primary vertex = ({},{},{})cm \t b field = {} tesla", pVtxPos.x, pVtxPos.y, pVtxPos.z, m_cfg.b_field/dd4hep::tesla); - - std::vector pi_index; - std::vector k_index; - std::vector p_index; + + std::vector hVec; hVec.clear(); + std::vector indexVec; indexVec.clear(); for (unsigned int i = 0; const auto& p : *rcparts) { - const auto pdg = p.getPDG(); - if(abs(pdg) == 211) pi_index.push_back(i); - if(abs(pdg) == 321) k_index.push_back(i); - if(abs(pdg) == 2212) p_index.push_back(i); - ++i; + if(p.getCharge()==0) continue; + Helix h(p, m_cfg.b_field); + double dca = h.distance(pVtxPos) * edm4eic::unit::cm; + if( dca < m_cfg.minDca ) continue; + + hVec.push_back(h); + indexVec.push_back(i); + ++i; } - info("\t Array sizes: pions = {}, kaons = {}, protons = {}", pi_index.size(), k_index.size(), p_index.size()); + if( hVec.size() != indexVec.size() ) return; + + debug("\t Vector size {}, {}", hVec.size(), indexVec.size()); + + for (unsigned int i1 = 0; i1 < hVec.size(); ++i1) { + for (unsigned int i2 = i1 + 1; i2 < hVec.size(); ++i2) { + const auto& p1 = (*rcparts)[indexVec[i1]]; + const auto& p2 = (*rcparts)[indexVec[i2]]; - for (unsigned int i1 = 0; i1 < pi_index.size(); ++i1) { - for (unsigned int i2 = i1 + 1; i2 < pi_index.size(); ++i2) { - const auto& p1 = (*rcparts)[i1]; - const auto& p2 = (*rcparts)[i2]; - - if (p1.getCharge() + p2.getCharge() != 0) continue; + if ( !(m_cfg.unlikesign && p1.getCharge() + p2.getCharge() == 0) ) continue; - Helix h1obj(p1, m_cfg.b_field); Helix& h1 = h1obj; - Helix h2obj(p2, m_cfg.b_field); Helix& h2 = h2obj; + const auto& h1 = hVec[i1]; + const auto& h2 = hVec[i2]; // Helix function uses cm unit double dca1 = h1.distance(pVtxPos) * edm4eic::unit::cm; double dca2 = h2.distance(pVtxPos) * edm4eic::unit::cm; - debug("\t dca1 = {}, dca2 = {}", dca1, dca2); - if( dca1 < m_cfg.minDca1 || dca2 < m_cfg.minDca2 ) continue; - + if( dca1 < m_cfg.minDca || dca2 < m_cfg.minDca ) continue; + std::pair const ss = h1.pathLengths(h2); edm4hep::Vector3f h1AtDcaTo2 = h1.at(ss.first); edm4hep::Vector3f h2AtDcaTo1 = h2.at(ss.second); double dca12 = edm4hep::utils::magnitude(h1AtDcaTo2 - h2AtDcaTo1) * edm4eic::unit::cm; + if(std::isnan(dca12)) continue; if( dca12 > m_cfg.maxDca12 ) continue; edm4hep::Vector3f pairPos = 0.5*(h1AtDcaTo2 + h2AtDcaTo1); @@ -90,8 +97,8 @@ void SecondaryVerticesHelix::process(const SecondaryVerticesHelix::Input& input, edm4hep::Vector3f h2MomAtDca = h2.momentumAt(ss.second, m_cfg.b_field); edm4hep::Vector3f pairMom = h1MomAtDca + h2MomAtDca; - double e1 = std::hypot(edm4hep::utils::magnitude(h1MomAtDca), particleSvc.particle(211).mass); - double e2 = std::hypot(edm4hep::utils::magnitude(h2MomAtDca), particleSvc.particle(211).mass); + double e1 = std::hypot(edm4hep::utils::magnitude(h1MomAtDca), particleSvc.particle(p1.getPDG()).mass); + double e2 = std::hypot(edm4hep::utils::magnitude(h2MomAtDca), particleSvc.particle(p2.getPDG()).mass); double pairE = e1+e2; double pairP = edm4hep::utils::magnitude(pairMom); @@ -102,6 +109,11 @@ void SecondaryVerticesHelix::process(const SecondaryVerticesHelix::Input& input, double beta = edm4hep::utils::magnitude(pairMom)/pairE; double time = edm4hep::utils::magnitude(pairPos - pVtxPos)/(beta*dd4hep::c_light); + edm4hep::Vector3f dL = pairPos - pVtxPos; // in cm + edm4hep::Vector3f decayL(dL.x*edm4eic::unit::cm, dL.y*edm4eic::unit::cm, dL.z*edm4eic::unit::cm); + double dca2pv = edm4hep::utils::magnitude(decayL) * sin(angle); + if(dca2pv > m_cfg.maxDca) continue; + auto v0 = out_secondary_vertices->create(); v0.setType(2); // 2 for secondary v0.setPosition({(float)(pairPos.x * edm4eic::unit::cm / edm4eic::unit::mm), @@ -109,7 +121,7 @@ void SecondaryVerticesHelix::process(const SecondaryVerticesHelix::Input& input, (float)(pairPos.z * edm4eic::unit::cm / edm4eic::unit::mm), (float)time}); v0.addToAssociatedParticles(p1); - v0.addToAssociatedParticles(p2); + v0.addToAssociatedParticles(p2); info("One secondary vertex found at (x,y,z) = ({}, {}, {}) mm.", pairPos.x * edm4eic::unit::cm / edm4eic::unit::mm, diff --git a/src/algorithms/reco/SecondaryVerticesHelix.h b/src/algorithms/reco/SecondaryVerticesHelix.h index 378cb45558..dc3ece37a2 100644 --- a/src/algorithms/reco/SecondaryVerticesHelix.h +++ b/src/algorithms/reco/SecondaryVerticesHelix.h @@ -1,5 +1,5 @@ // SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2024 Daniel Brandenburg, Xin Dong +// Copyright (C) 2025 Xin Dong #pragma once diff --git a/src/algorithms/reco/SecondaryVerticesHelixConfig.h b/src/algorithms/reco/SecondaryVerticesHelixConfig.h index 00e47c6ae1..3321e9fab5 100644 --- a/src/algorithms/reco/SecondaryVerticesHelixConfig.h +++ b/src/algorithms/reco/SecondaryVerticesHelixConfig.h @@ -12,11 +12,11 @@ namespace eicrecon { struct SecondaryVerticesHelixConfig { float b_field = -1.7*dd4hep::tesla; - float minDca1 = 0.03*edm4eic::unit::mm; // mm, daughter1 to pVtx - float minDca2 = 0.03*edm4eic::unit::mm; // mm, daughter2 to pVtx - float maxDca12 = 0.3*edm4eic::unit::mm; // mm, dca between daughter 1 and 2 - float maxDca = 0.3*edm4eic::unit::mm; // mm, dca of V0 to pVtx - float minCostheta = 0.95; // costheta, theta: angle of V0 decay direction and momentum + bool unlikesign = true; + float minDca = 0.03*edm4eic::unit::mm; // mm, daughter to pVtx + float maxDca12 = 1.*edm4eic::unit::mm; // mm, dca between daughter 1 and 2 + float maxDca = 1.*edm4eic::unit::mm; // mm, dca of V0 to pVtx + float minCostheta = 0.8; // costheta, theta: angle of V0 decay direction and momentum }; } // namespace eicrecon From 21912b102a1410c24de7747fe2c8d216c85baa4f Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 20 Oct 2025 21:15:01 +0000 Subject: [PATCH 07/17] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/algorithms/reco/Helix.cc | 1029 ++++++++--------- src/algorithms/reco/Helix.h | 338 +++--- src/algorithms/reco/SecondaryVerticesHelix.cc | 160 +-- src/algorithms/reco/SecondaryVerticesHelix.h | 11 +- .../reco/SecondaryVerticesHelixConfig.h | 10 +- .../reco/SecondaryVerticesHelix_factory.h | 3 +- src/global/reco/reco.cc | 4 +- 7 files changed, 760 insertions(+), 795 deletions(-) diff --git a/src/algorithms/reco/Helix.cc b/src/algorithms/reco/Helix.cc index eed158857e..6664fb3f8e 100644 --- a/src/algorithms/reco/Helix.cc +++ b/src/algorithms/reco/Helix.cc @@ -25,15 +25,13 @@ namespace eicrecon { const double Helix::NoSolution = 3.e+33; // basic constructor -Helix::Helix(double c, double d, double phase, const edm4hep::Vector3f& o, int h) -{ - setParameters(c, d, phase, o, h); +Helix::Helix(double c, double d, double phase, const edm4hep::Vector3f& o, int h) { + setParameters(c, d, phase, o, h); } // momentum in GeV, position in cm, B in Tesla -Helix::Helix(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q) -{ - setParameters(p, o, B, q); +Helix::Helix(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q) { + setParameters(p, o, B, q); } // construct using ReconstructParticle @@ -50,508 +48,481 @@ Helix::Helix(const edm4eic::ReconstructedParticle& p, const double b_field) { // construct using TrackParameteters void Helix::setParameters(const edm4eic::TrackParameters& trk, const double b_field) { - const auto mom = edm4hep::utils::sphericalToVector(1.0 / std::abs(trk.getQOverP()), trk.getTheta(), trk.getPhi()); + const auto mom = edm4hep::utils::sphericalToVector(1.0 / std::abs(trk.getQOverP()), + trk.getTheta(), trk.getPhi()); const auto charge = std::copysign(1., trk.getQOverP()); - const auto phi = trk.getPhi(); - const auto loc = trk.getLoc(); - edm4hep::Vector3f pos( -1. * loc.a * sin(phi), loc.a * cos(phi), loc.b); // PCA point + const auto phi = trk.getPhi(); + const auto loc = trk.getLoc(); + edm4hep::Vector3f pos(-1. * loc.a * sin(phi), loc.a * cos(phi), loc.b); // PCA point - setParameters(mom*edm4eic::unit::GeV/dd4hep::GeV, pos*edm4eic::unit::mm/edm4eic::unit::cm, b_field, charge); + setParameters(mom * edm4eic::unit::GeV / dd4hep::GeV, pos * edm4eic::unit::mm / edm4eic::unit::cm, + b_field, charge); } -void Helix::setParameters(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q) -{ - mH = (q*B <= 0) ? 1 : -1; - if(p.y == 0 && p.x == 0) - setPhase((M_PI/4)*(1-2.*mH)); - else - setPhase(atan2(p.y,p.x)-mH*M_PI/2); - setDipAngle(atan2(p.z,edm4hep::utils::magnitudeTransverse(p))); - mOrigin = o; - - double curvature_val = fabs((dd4hep::c_light*dd4hep::nanosecond/dd4hep::meter*q*B/dd4hep::tesla)/ - (edm4hep::utils::magnitude(p)/dd4hep::GeV*mCosDipAngle)/dd4hep::meter); - setCurvature(curvature_val); +void Helix::setParameters(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, + const int q) { + mH = (q * B <= 0) ? 1 : -1; + if (p.y == 0 && p.x == 0) + setPhase((M_PI / 4) * (1 - 2. * mH)); + else + setPhase(atan2(p.y, p.x) - mH * M_PI / 2); + setDipAngle(atan2(p.z, edm4hep::utils::magnitudeTransverse(p))); + mOrigin = o; + + double curvature_val = + fabs((dd4hep::c_light * dd4hep::nanosecond / dd4hep::meter * q * B / dd4hep::tesla) / + (edm4hep::utils::magnitude(p) / dd4hep::GeV * mCosDipAngle) / dd4hep::meter); + setCurvature(curvature_val); } -void Helix::setParameters(double c, double dip, double phase, - const edm4hep::Vector3f& o, int h) -{ - // - // The order in which the parameters are set is important - // since setCurvature might have to adjust the others. - // - mH = (h>=0) ? 1 : -1; // Default is: positive particle - // positive field - mOrigin = o; - setDipAngle(dip); - setPhase(phase); - - // - // Check for singularity and correct for negative curvature. - // May change mH and mPhase. Must therefore be set last. - // - setCurvature(c); - - // - // For the case B=0, h is ill defined. In the following we - // always assume h = +1. Since phase = psi - h * pi/2 - // we have to correct the phase in case h = -1. - // This assumes that the user uses the same h for phase - // as the one he passed to the constructor. - // - if (mSingularity && mH == -1) { - mH = +1; - setPhase(mPhase-M_PI); - } +void Helix::setParameters(double c, double dip, double phase, const edm4hep::Vector3f& o, int h) { + // + // The order in which the parameters are set is important + // since setCurvature might have to adjust the others. + // + mH = (h >= 0) ? 1 : -1; // Default is: positive particle + // positive field + mOrigin = o; + setDipAngle(dip); + setPhase(phase); + + // + // Check for singularity and correct for negative curvature. + // May change mH and mPhase. Must therefore be set last. + // + setCurvature(c); + + // + // For the case B=0, h is ill defined. In the following we + // always assume h = +1. Since phase = psi - h * pi/2 + // we have to correct the phase in case h = -1. + // This assumes that the user uses the same h for phase + // as the one he passed to the constructor. + // + if (mSingularity && mH == -1) { + mH = +1; + setPhase(mPhase - M_PI); + } } -void Helix::setCurvature(double val) -{ - if (val < 0) { - mCurvature = -val; - mH = -mH; - setPhase(mPhase+M_PI); - } - else - mCurvature = val; - - if (fabs(mCurvature) <= std::numeric_limits::epsilon()) - mSingularity = true; // straight line - else - mSingularity = false; // curved +void Helix::setCurvature(double val) { + if (val < 0) { + mCurvature = -val; + mH = -mH; + setPhase(mPhase + M_PI); + } else + mCurvature = val; + + if (fabs(mCurvature) <= std::numeric_limits::epsilon()) + mSingularity = true; // straight line + else + mSingularity = false; // curved } -void Helix::setPhase(double val) -{ - mPhase = val; - mCosPhase = cos(mPhase); - mSinPhase = sin(mPhase); - if (fabs(mPhase) > M_PI) - mPhase = atan2(mSinPhase, mCosPhase); // force range [-pi,pi] +void Helix::setPhase(double val) { + mPhase = val; + mCosPhase = cos(mPhase); + mSinPhase = sin(mPhase); + if (fabs(mPhase) > M_PI) + mPhase = atan2(mSinPhase, mCosPhase); // force range [-pi,pi] } -void Helix::setDipAngle(double val) -{ - mDipAngle = val; - mCosDipAngle = cos(mDipAngle); - mSinDipAngle = sin(mDipAngle); +void Helix::setDipAngle(double val) { + mDipAngle = val; + mCosDipAngle = cos(mDipAngle); + mSinDipAngle = sin(mDipAngle); } -double Helix::xcenter() const -{ - if (mSingularity) - return 0; - else - return mOrigin.x-mCosPhase/mCurvature; +double Helix::xcenter() const { + if (mSingularity) + return 0; + else + return mOrigin.x - mCosPhase / mCurvature; } -double Helix::ycenter() const -{ - if (mSingularity) - return 0; - else - return mOrigin.y-mSinPhase/mCurvature; +double Helix::ycenter() const { + if (mSingularity) + return 0; + else + return mOrigin.y - mSinPhase / mCurvature; } -double Helix::fudgePathLength(const edm4hep::Vector3f& p) const -{ - double s; - double dx = p.x-mOrigin.x; - double dy = p.y-mOrigin.y; - - if (mSingularity) { - s = (dy*mCosPhase - dx*mSinPhase)/mCosDipAngle; - } - else { - s = atan2(dy*mCosPhase - dx*mSinPhase, - 1/mCurvature + dx*mCosPhase+dy*mSinPhase)/ - (mH*mCurvature*mCosDipAngle); - } - return s; +double Helix::fudgePathLength(const edm4hep::Vector3f& p) const { + double s; + double dx = p.x - mOrigin.x; + double dy = p.y - mOrigin.y; + + if (mSingularity) { + s = (dy * mCosPhase - dx * mSinPhase) / mCosDipAngle; + } else { + s = atan2(dy * mCosPhase - dx * mSinPhase, 1 / mCurvature + dx * mCosPhase + dy * mSinPhase) / + (mH * mCurvature * mCosDipAngle); + } + return s; } -double Helix::distance(const edm4hep::Vector3f& p, bool scanPeriods) const -{ - return edm4hep::utils::magnitude(this->at(pathLength(p,scanPeriods))-p); +double Helix::distance(const edm4hep::Vector3f& p, bool scanPeriods) const { + return edm4hep::utils::magnitude(this->at(pathLength(p, scanPeriods)) - p); } -double Helix::pathLength(const edm4hep::Vector3f& p, bool scanPeriods) const -{ +double Helix::pathLength(const edm4hep::Vector3f& p, bool scanPeriods) const { + // + // Returns the path length at the distance of closest + // approach between the helix and point p. + // For the case of B=0 (straight line) the path length + // can be calculated analytically. For B>0 there is + // unfortunately no easy solution to the problem. + // Here we use the Newton method to find the root of the + // referring equation. The 'fudgePathLength' serves + // as a starting value. + // + double s; + double dx = p.x - mOrigin.x; + double dy = p.y - mOrigin.y; + double dz = p.z - mOrigin.z; + + if (mSingularity) { + s = mCosDipAngle * (mCosPhase * dy - mSinPhase * dx) + mSinDipAngle * dz; + } else { // + const double MaxPrecisionNeeded = edm4eic::unit::um; + const int MaxIterations = 100; + // - // Returns the path length at the distance of closest - // approach between the helix and point p. - // For the case of B=0 (straight line) the path length - // can be calculated analytically. For B>0 there is - // unfortunately no easy solution to the problem. - // Here we use the Newton method to find the root of the - // referring equation. The 'fudgePathLength' serves - // as a starting value. + // The math is taken from Maple with C(expr,optimized) and + // some hand-editing. It is not very nice but efficient. // - double s; - double dx = p.x-mOrigin.x; - double dy = p.y-mOrigin.y; - double dz = p.z-mOrigin.z; + double t34 = mCurvature * mCosDipAngle * mCosDipAngle; + double t41 = mSinDipAngle * mSinDipAngle; + double t6, t7, t11, t12, t19; - if (mSingularity) { - s = mCosDipAngle*(mCosPhase*dy-mSinPhase*dx) + - mSinDipAngle*dz; + // + // Get a first guess by using the dca in 2D. Since + // in some extreme cases we might be off by n periods + // we add (subtract) periods in case we get any closer. + // + s = fudgePathLength(p); + + if (scanPeriods) { + double ds = period(); + int j, jmin = 0; + double d, dmin = edm4hep::utils::magnitude(at(s) - p); + for (j = 1; j < MaxIterations; j++) { + if ((d = edm4hep::utils::magnitude(at(s + j * ds) - p)) < dmin) { + dmin = d; + jmin = j; + } else + break; + } + for (j = -1; -j < MaxIterations; j--) { + if ((d = edm4hep::utils::magnitude(at(s + j * ds) - p)) < dmin) { + dmin = d; + jmin = j; + } else + break; + } + if (jmin) + s += jmin * ds; } - else { // - const double MaxPrecisionNeeded = edm4eic::unit::um; - const int MaxIterations = 100; - - // - // The math is taken from Maple with C(expr,optimized) and - // some hand-editing. It is not very nice but efficient. - // - double t34 = mCurvature*mCosDipAngle*mCosDipAngle; - double t41 = mSinDipAngle*mSinDipAngle; - double t6, t7, t11, t12, t19; - - // - // Get a first guess by using the dca in 2D. Since - // in some extreme cases we might be off by n periods - // we add (subtract) periods in case we get any closer. - // - s = fudgePathLength(p); - - if (scanPeriods) { - double ds = period(); - int j, jmin = 0; - double d, dmin = edm4hep::utils::magnitude(at(s) - p); - for(j=1; j::max(); - else - return fabs(2*M_PI/(mH*mCurvature*mCosDipAngle)); +double Helix::period() const { + if (mSingularity) + return std::numeric_limits::max(); + else + return fabs(2 * M_PI / (mH * mCurvature * mCosDipAngle)); } -std::pair Helix::pathLength(double r) const -{ - std::pair value; - std::pair VALUE(999999999.,999999999.); +std::pair Helix::pathLength(double r) const { + std::pair value; + std::pair VALUE(999999999., 999999999.); + // + // The math is taken from Maple with C(expr,optimized) and + // some hand-editing. It is not very nice but efficient. + // 'first' is the smallest of the two solutions (may be negative) + // 'second' is the other. + // + if (mSingularity) { + double t1 = mCosDipAngle * (mOrigin.x * mSinPhase - mOrigin.y * mCosPhase); + double t12 = mOrigin.y * mOrigin.y; + double t13 = mCosPhase * mCosPhase; + double t15 = r * r; + double t16 = mOrigin.x * mOrigin.x; + double t20 = + -mCosDipAngle * mCosDipAngle * + (2.0 * mOrigin.x * mSinPhase * mOrigin.y * mCosPhase + t12 - t12 * t13 - t15 + t13 * t16); + if (t20 < 0.) + return VALUE; + t20 = ::sqrt(t20); + value.first = (t1 - t20) / (mCosDipAngle * mCosDipAngle); + value.second = (t1 + t20) / (mCosDipAngle * mCosDipAngle); + } else { + double t1 = mOrigin.y * mCurvature; + double t2 = mSinPhase; + double t3 = mCurvature * mCurvature; + double t4 = mOrigin.y * t2; + double t5 = mCosPhase; + double t6 = mOrigin.x * t5; + double t8 = mOrigin.x * mOrigin.x; + double t11 = mOrigin.y * mOrigin.y; + double t14 = r * r; + double t15 = t14 * mCurvature; + double t17 = t8 * t8; + double t19 = t11 * t11; + double t21 = t11 * t3; + double t23 = t5 * t5; + double t32 = t14 * t14; + double t35 = t14 * t3; + double t38 = 8.0 * t4 * t6 - 4.0 * t1 * t2 * t8 - 4.0 * t11 * mCurvature * t6 + 4.0 * t15 * t6 + + t17 * t3 + t19 * t3 + 2.0 * t21 * t8 + 4.0 * t8 * t23 - + 4.0 * t8 * mOrigin.x * mCurvature * t5 - 4.0 * t11 * t23 - + 4.0 * t11 * mOrigin.y * mCurvature * t2 + 4.0 * t11 - 4.0 * t14 + t32 * t3 + + 4.0 * t15 * t4 - 2.0 * t35 * t11 - 2.0 * t35 * t8; + double t40 = (-t3 * t38); + if (t40 < 0.) + return VALUE; + t40 = ::sqrt(t40); + + double t43 = mOrigin.x * mCurvature; + double t45 = 2.0 * t5 - t35 + t21 + 2.0 - 2.0 * t1 * t2 - 2.0 * t43 - 2.0 * t43 * t5 + t8 * t3; + double t46 = mH * mCosDipAngle * mCurvature; + + value.first = (-mPhase + 2.0 * atan((-2.0 * t1 + 2.0 * t2 + t40) / t45)) / t46; + value.second = -(mPhase + 2.0 * atan((2.0 * t1 - 2.0 * t2 + t40) / t45)) / t46; + // - // The math is taken from Maple with C(expr,optimized) and - // some hand-editing. It is not very nice but efficient. - // 'first' is the smallest of the two solutions (may be negative) - // 'second' is the other. + // Solution can be off by +/- one period, select smallest // - if (mSingularity) { - double t1 = mCosDipAngle*(mOrigin.x*mSinPhase-mOrigin.y*mCosPhase); - double t12 = mOrigin.y*mOrigin.y; - double t13 = mCosPhase*mCosPhase; - double t15 = r*r; - double t16 = mOrigin.x*mOrigin.x; - double t20 = -mCosDipAngle*mCosDipAngle*(2.0*mOrigin.x*mSinPhase*mOrigin.y*mCosPhase + - t12-t12*t13-t15+t13*t16); - if (t20<0.) return VALUE; - t20 = ::sqrt(t20); - value.first = (t1-t20)/(mCosDipAngle*mCosDipAngle); - value.second = (t1+t20)/(mCosDipAngle*mCosDipAngle); + double p = period(); + if (!std::isnan(value.first)) { + if (fabs(value.first - p) < fabs(value.first)) + value.first = value.first - p; + else if (fabs(value.first + p) < fabs(value.first)) + value.first = value.first + p; } - else { - double t1 = mOrigin.y*mCurvature; - double t2 = mSinPhase; - double t3 = mCurvature*mCurvature; - double t4 = mOrigin.y*t2; - double t5 = mCosPhase; - double t6 = mOrigin.x*t5; - double t8 = mOrigin.x*mOrigin.x; - double t11 = mOrigin.y*mOrigin.y; - double t14 = r*r; - double t15 = t14*mCurvature; - double t17 = t8*t8; - double t19 = t11*t11; - double t21 = t11*t3; - double t23 = t5*t5; - double t32 = t14*t14; - double t35 = t14*t3; - double t38 = 8.0*t4*t6 - 4.0*t1*t2*t8 - 4.0*t11*mCurvature*t6 + - 4.0*t15*t6 + t17*t3 + t19*t3 + 2.0*t21*t8 + 4.0*t8*t23 - - 4.0*t8*mOrigin.x*mCurvature*t5 - 4.0*t11*t23 - - 4.0*t11*mOrigin.y*mCurvature*t2 + 4.0*t11 - 4.0*t14 + - t32*t3 + 4.0*t15*t4 - 2.0*t35*t11 - 2.0*t35*t8; - double t40 = (-t3*t38); - if (t40<0.) return VALUE; - t40 = ::sqrt(t40); - - double t43 = mOrigin.x*mCurvature; - double t45 = 2.0*t5 - t35 + t21 + 2.0 - 2.0*t1*t2 -2.0*t43 - 2.0*t43*t5 + t8*t3; - double t46 = mH*mCosDipAngle*mCurvature; - - value.first = (-mPhase + 2.0*atan((-2.0*t1 + 2.0*t2 + t40)/t45))/t46; - value.second = -(mPhase + 2.0*atan((2.0*t1 - 2.0*t2 + t40)/t45))/t46; - - // - // Solution can be off by +/- one period, select smallest - // - double p = period(); - if (! std::isnan(value.first)) { - if (fabs(value.first-p) < fabs(value.first)) value.first = value.first-p; - else if (fabs(value.first+p) < fabs(value.first)) value.first = value.first+p; - } - if (! std::isnan(value.second)) { - if (fabs(value.second-p) < fabs(value.second)) value.second = value.second-p; - else if (fabs(value.second+p) < fabs(value.second)) value.second = value.second+p; - } + if (!std::isnan(value.second)) { + if (fabs(value.second - p) < fabs(value.second)) + value.second = value.second - p; + else if (fabs(value.second + p) < fabs(value.second)) + value.second = value.second + p; } - if (value.first > value.second) - std::swap(value.first,value.second); - return(value); + } + if (value.first > value.second) + std::swap(value.first, value.second); + return (value); } -std::pair Helix::pathLength(double r, double x, double y) -{ - double x0 = mOrigin.x; - double y0 = mOrigin.y; - mOrigin.x = x0-x; - mOrigin.y = y0-y; - std::pair result = this->pathLength(r); - mOrigin.x = x0; - mOrigin.y = y0; - return result; +std::pair Helix::pathLength(double r, double x, double y) { + double x0 = mOrigin.x; + double y0 = mOrigin.y; + mOrigin.x = x0 - x; + mOrigin.y = y0 - y; + std::pair result = this->pathLength(r); + mOrigin.x = x0; + mOrigin.y = y0; + return result; } -double Helix::pathLength(const edm4hep::Vector3f& r, - const edm4hep::Vector3f& n) const -{ +double Helix::pathLength(const edm4hep::Vector3f& r, const edm4hep::Vector3f& n) const { + // + // Vector 'r' defines the position of the center and + // vector 'n' the normal vector of the plane. + // For a straight line there is a simple analytical + // solution. For curvatures > 0 the root is determined + // by Newton method. In case no valid s can be found + // the max. largest value for s is returned. + // + double s; + + if (mSingularity) { + double t = n.z * mSinDipAngle + n.y * mCosDipAngle * mCosPhase - n.x * mCosDipAngle * mSinPhase; + if (t == 0) + s = NoSolution; + else + s = ((r - mOrigin) * n) / t; + } else { + const double MaxPrecisionNeeded = edm4eic::unit::um; + const int MaxIterations = 20; + + double A = mCurvature * ((mOrigin - r) * n) - n.x * mCosPhase - n.y * mSinPhase; + double t = mH * mCurvature * mCosDipAngle; + double u = n.z * mCurvature * mSinDipAngle; + + double a, f, fp; + double sOld = s = 0; + double shiftOld = 0; + double shift; + // (cos(angMax)-1)/angMax = 0.1 + const double angMax = 0.21; + double deltas = fabs(angMax / (mCurvature * mCosDipAngle)); + // dampingFactor = exp(-0.5); + // double dampingFactor = 0.60653; + int i; + + for (i = 0; i < MaxIterations; i++) { + a = t * s + mPhase; + double sina = sin(a); + double cosa = cos(a); + f = A + n.x * cosa + n.y * sina + u * s; + fp = -n.x * sina * t + n.y * cosa * t + u; + if (fabs(fp) * deltas <= fabs(f)) { //too big step + int sgn = 1; + if (fp < 0.) + sgn = -sgn; + if (f < 0.) + sgn = -sgn; + shift = sgn * deltas; + if (shift < 0) + shift *= 0.9; // don't get stuck shifting +/-deltas + } else { + shift = f / fp; + } + s -= shift; + shiftOld = shift; + if (fabs(sOld - s) < MaxPrecisionNeeded) + break; + sOld = s; + } + if (i == MaxIterations) + return NoSolution; + } + return s; +} + +std::pair Helix::pathLengths(const Helix& h, double minStepSize, + double minRange) const { + // + // Cannot handle case where one is a helix + // and the other one is a straight line. + // + if (mSingularity != h.mSingularity) + return std::pair(NoSolution, NoSolution); + + double s1, s2; + + if (mSingularity) { // - // Vector 'r' defines the position of the center and - // vector 'n' the normal vector of the plane. - // For a straight line there is a simple analytical - // solution. For curvatures > 0 the root is determined - // by Newton method. In case no valid s can be found - // the max. largest value for s is returned. + // Analytic solution // - double s; + edm4hep::Vector3f dv = h.mOrigin - mOrigin; + edm4hep::Vector3f a(-mCosDipAngle * mSinPhase, mCosDipAngle * mCosPhase, mSinDipAngle); + edm4hep::Vector3f b(-h.mCosDipAngle * h.mSinPhase, h.mCosDipAngle * h.mCosPhase, + h.mSinDipAngle); + double ab = a * b; + double g = dv * a; + double k = dv * b; + s2 = (k - ab * g) / (ab * ab - 1.); + s1 = g + s2 * ab; + return std::pair(s1, s2); + } else { + // + // First step: get dca in the xy-plane as start value + // + double dx = h.xcenter() - xcenter(); + double dy = h.ycenter() - ycenter(); + double dd = ::sqrt(dx * dx + dy * dy); + double r1 = 1 / curvature(); + double r2 = 1 / h.curvature(); - if (mSingularity) { - double t = n.z*mSinDipAngle + - n.y*mCosDipAngle*mCosPhase - - n.x*mCosDipAngle*mSinPhase; - if (t == 0) - s = NoSolution; - else - s = ((r - mOrigin)*n)/t; - } - else { - const double MaxPrecisionNeeded = edm4eic::unit::um; - const int MaxIterations = 20; - - double A = mCurvature*((mOrigin - r)*n) - - n.x*mCosPhase - - n.y*mSinPhase; - double t = mH*mCurvature*mCosDipAngle; - double u = n.z*mCurvature*mSinDipAngle; - - double a, f, fp; - double sOld = s = 0; - double shiftOld = 0; - double shift; -// (cos(angMax)-1)/angMax = 0.1 - const double angMax = 0.21; - double deltas = fabs(angMax/(mCurvature*mCosDipAngle)); -// dampingFactor = exp(-0.5); -// double dampingFactor = 0.60653; - int i; - - for (i=0; i dd ? -1 : 1); // set -1 when *this* helix is + // completely contained in the other + x = xcenter() + rsign * r1 * dx / dd; + y = ycenter() + rsign * r1 * dy / dd; + s = pathLength(x, y); } - return s; -} -std::pair -Helix::pathLengths(const Helix& h, double minStepSize, double minRange) const -{ // - // Cannot handle case where one is a helix - // and the other one is a straight line. + // Second step: scan in decreasing intervals around seed 's' + // minRange and minStepSize are passed as arguments to the method. + // They have default values defined in the header file. // - if (mSingularity != h.mSingularity) - return std::pair(NoSolution, NoSolution); - - double s1, s2; - - if (mSingularity) { - // - // Analytic solution - // - edm4hep::Vector3f dv = h.mOrigin - mOrigin; - edm4hep::Vector3f a(-mCosDipAngle*mSinPhase, - mCosDipAngle*mCosPhase, - mSinDipAngle); - edm4hep::Vector3f b(-h.mCosDipAngle*h.mSinPhase, - h.mCosDipAngle*h.mCosPhase, - h.mSinDipAngle); - double ab = a*b; - double g = dv*a; - double k = dv*b; - s2 = (k-ab*g)/(ab*ab-1.); - s1 = g+s2*ab; - return std::pair(s1, s2); - } - else { - // - // First step: get dca in the xy-plane as start value - // - double dx = h.xcenter() - xcenter(); - double dy = h.ycenter() - ycenter(); - double dd = ::sqrt(dx*dx + dy*dy); - double r1 = 1/curvature(); - double r2 = 1/h.curvature(); - - double cosAlpha = (r1*r1 + dd*dd - r2*r2)/(2*r1*dd); - - double s; - double x, y; - if (fabs(cosAlpha) < 1) { // two solutions - double sinAlpha = sin(acos(cosAlpha)); - x = xcenter() + r1*(cosAlpha*dx - sinAlpha*dy)/dd; - y = ycenter() + r1*(sinAlpha*dx + cosAlpha*dy)/dd; - s = pathLength(x, y); - x = xcenter() + r1*(cosAlpha*dx + sinAlpha*dy)/dd; - y = ycenter() + r1*(cosAlpha*dy - sinAlpha*dx)/dd; - double a = pathLength(x, y); - if (h.distance(at(a)) < h.distance(at(s))) s = a; - } - else { // no intersection (or exactly one) - int rsign = ((r2-r1) > dd ? -1 : 1); // set -1 when *this* helix is - // completely contained in the other - x = xcenter() + rsign*r1*dx/dd; - y = ycenter() + rsign*r1*dy/dd; - s = pathLength(x, y); - } - - // - // Second step: scan in decreasing intervals around seed 's' - // minRange and minStepSize are passed as arguments to the method. - // They have default values defined in the header file. - // - double dmin = h.distance(at(s)); - double range = std::max(2*dmin, minRange); - double ds = range/10; - double slast=-999999, ss, d; - s1 = s - range/2.; - s2 = s + range/2.; - - while (ds > minStepSize) { - for (ss=s1; ss minStepSize) { + for (ss = s1; ss < s2 + ds; ss += ds) { + d = h.distance(at(ss)); + if (d < dmin) { + dmin = d; + s = ss; } - return std::pair(s, h.pathLength(at(s))); + slast = ss; + } + // + // In the rare cases where the minimum is at the + // the border of the current range we shift the range + // and start all over, i.e we do not decrease 'ds'. + // Else we decrease the search intervall around the + // current minimum and redo the scan in smaller steps. + // + if (s == s1) { + d = 0.8 * (s2 - s1); + s1 -= d; + s2 -= d; + } else if (s == slast) { + d = 0.8 * (s2 - s1); + s1 += d; + s2 += d; + } else { + s1 = s - ds; + s2 = s + ds; + ds /= 10; + } } + return std::pair(s, h.pathLength(at(s))); + } } - -void Helix::moveOrigin(double s) -{ - if (mSingularity) - mOrigin = at(s); - else { - edm4hep::Vector3f newOrigin = at(s); - double newPhase = atan2(newOrigin.y - ycenter(), - newOrigin.x - xcenter()); - mOrigin = newOrigin; - setPhase(newPhase); - } +void Helix::moveOrigin(double s) { + if (mSingularity) + mOrigin = at(s); + else { + edm4hep::Vector3f newOrigin = at(s); + double newPhase = atan2(newOrigin.y - ycenter(), newOrigin.x - xcenter()); + mOrigin = newOrigin; + setPhase(newPhase); + } } /* int operator== (const Helix& a, const Helix& b) @@ -560,100 +531,88 @@ int operator== (const Helix& a, const Helix& b) // Checks for numerical identity only ! // return (a.origin() == b.origin() && - a.dipAngle() == b.dipAngle() && - a.curvature() == b.curvature() && - a.phase() == b.phase() && - a.h() == b.h()); + a.dipAngle() == b.dipAngle() && + a.curvature() == b.curvature() && + a.phase() == b.phase() && + a.h() == b.h()); } int operator!= (const Helix& a, const Helix& b) {return !(a == b);} */ //std::ostream& operator<<(std::ostream& os, const Helix& h) -void Helix::Print() const -{ - std::cout << "(" - << "curvature = " << mCurvature << ", " - << "dip angle = " << mDipAngle << ", " - << "phase = " << mPhase << ", " - << "h = " << mH << ", " - << "origin = " << mOrigin.x << " " << mOrigin.y << " " << mOrigin.z << ")" - << std::endl; +void Helix::Print() const { + std::cout << "(" + << "curvature = " << mCurvature << ", " + << "dip angle = " << mDipAngle << ", " + << "phase = " << mPhase << ", " + << "h = " << mH << ", " + << "origin = " << mOrigin.x << " " << mOrigin.y << " " << mOrigin.z << ")" << std::endl; } -edm4hep::Vector3f Helix::momentum(double B) const -{ - if (mSingularity) - return(edm4hep::Vector3f(0,0,0)); - else { - double pt = edm4eic::unit::GeV*fabs(dd4hep::c_light*dd4hep::nanosecond/dd4hep::meter*B/dd4hep::tesla)/(fabs(mCurvature)*dd4hep::meter); - - return (edm4hep::Vector3f(pt*cos(mPhase+mH*M_PI/2), // pos part pos field - pt*sin(mPhase+mH*M_PI/2), - pt*tan(mDipAngle))); - } -} +edm4hep::Vector3f Helix::momentum(double B) const { + if (mSingularity) + return (edm4hep::Vector3f(0, 0, 0)); + else { + double pt = edm4eic::unit::GeV * + fabs(dd4hep::c_light * dd4hep::nanosecond / dd4hep::meter * B / dd4hep::tesla) / + (fabs(mCurvature) * dd4hep::meter); -edm4hep::Vector3f Helix::momentumAt(double S, double B) const -{ - // Obtain phase-shifted momentum from phase-shift of origin - Helix tmp(*this); - tmp.moveOrigin(S); - return tmp.momentum(B); + return (edm4hep::Vector3f(pt * cos(mPhase + mH * M_PI / 2), // pos part pos field + pt * sin(mPhase + mH * M_PI / 2), pt * tan(mDipAngle))); + } } -int Helix::charge(double B) const -{ - return (B > 0 ? -mH : mH); +edm4hep::Vector3f Helix::momentumAt(double S, double B) const { + // Obtain phase-shifted momentum from phase-shift of origin + Helix tmp(*this); + tmp.moveOrigin(S); + return tmp.momentum(B); } -double Helix::geometricSignedDistance(double x, double y) -{ - // Geometric signed distance - double thePath = this->pathLength(x,y); - edm4hep::Vector3f DCA2dPosition = this->at(thePath); - DCA2dPosition.z = 0; - edm4hep::Vector3f position(x,y,0); - edm4hep::Vector3f DCAVec = (DCA2dPosition-position); - edm4hep::Vector3f momVec; - // Deal with straight tracks - if (this->mSingularity) { - momVec = this->at(1)- this->at(0); - momVec.z = 0; - } - else { - momVec = this->momentumAt(thePath,1./dd4hep::tesla); // Don't care about Bmag. Helicity is what matters. - momVec.z = 0; - } - - double cross = DCAVec.x*momVec.y - DCAVec.y*momVec.x; - double theSign = (cross>=0) ? 1. : -1.; - return theSign*edm4hep::utils::magnitudeTransverse(DCAVec); +int Helix::charge(double B) const { return (B > 0 ? -mH : mH); } + +double Helix::geometricSignedDistance(double x, double y) { + // Geometric signed distance + double thePath = this->pathLength(x, y); + edm4hep::Vector3f DCA2dPosition = this->at(thePath); + DCA2dPosition.z = 0; + edm4hep::Vector3f position(x, y, 0); + edm4hep::Vector3f DCAVec = (DCA2dPosition - position); + edm4hep::Vector3f momVec; + // Deal with straight tracks + if (this->mSingularity) { + momVec = this->at(1) - this->at(0); + momVec.z = 0; + } else { + momVec = this->momentumAt( + thePath, 1. / dd4hep::tesla); // Don't care about Bmag. Helicity is what matters. + momVec.z = 0; + } + + double cross = DCAVec.x * momVec.y - DCAVec.y * momVec.x; + double theSign = (cross >= 0) ? 1. : -1.; + return theSign * edm4hep::utils::magnitudeTransverse(DCAVec); } -double Helix::curvatureSignedDistance(double x, double y) -{ - // Protect against mH = 0 or zero field - if (this->mSingularity || abs(this->mH)<=0) { - return (this->geometricSignedDistance(x,y)); - } - else { - return (this->geometricSignedDistance(x,y))/(this->mH); - } - +double Helix::curvatureSignedDistance(double x, double y) { + // Protect against mH = 0 or zero field + if (this->mSingularity || abs(this->mH) <= 0) { + return (this->geometricSignedDistance(x, y)); + } else { + return (this->geometricSignedDistance(x, y)) / (this->mH); + } } -double Helix::geometricSignedDistance(const edm4hep::Vector3f& pos) -{ - double sdca2d = this->geometricSignedDistance(pos.x,pos.y); - double theSign = (sdca2d>=0) ? 1. : -1.; - return (this->distance(pos))*theSign; +double Helix::geometricSignedDistance(const edm4hep::Vector3f& pos) { + double sdca2d = this->geometricSignedDistance(pos.x, pos.y); + double theSign = (sdca2d >= 0) ? 1. : -1.; + return (this->distance(pos)) * theSign; } -double Helix::curvatureSignedDistance(const edm4hep::Vector3f& pos) -{ - double sdca2d = this->curvatureSignedDistance(pos.x,pos.y); - double theSign = (sdca2d>=0) ? 1. : -1.; - return (this->distance(pos))*theSign; +double Helix::curvatureSignedDistance(const edm4hep::Vector3f& pos) { + double sdca2d = this->curvatureSignedDistance(pos.x, pos.y); + double theSign = (sdca2d >= 0) ? 1. : -1.; + return (this->distance(pos)) * theSign; } -} // namespace eicrecon \ No newline at end of file +} // namespace eicrecon diff --git a/src/algorithms/reco/Helix.h b/src/algorithms/reco/Helix.h index 37a281d6b5..d68354b7b1 100644 --- a/src/algorithms/reco/Helix.h +++ b/src/algorithms/reco/Helix.h @@ -21,119 +21,118 @@ namespace eicrecon { class Helix { - bool mSingularity; // true for straight line case (B=0) - edm4hep::Vector3f mOrigin; - double mDipAngle; - double mCurvature; - double mPhase; - int mH; // -sign(q*B); - - double mCosDipAngle; - double mSinDipAngle; - double mCosPhase; - double mSinPhase; + bool mSingularity; // true for straight line case (B=0) + edm4hep::Vector3f mOrigin; + double mDipAngle; + double mCurvature; + double mPhase; + int mH; // -sign(q*B); + + double mCosDipAngle; + double mSinDipAngle; + double mCosPhase; + double mSinPhase; public: - /// curvature, dip angle, phase, origin, h - Helix(const double c, const double dip, const double phase, const edm4hep::Vector3f& o, const int h=-1); - - /// momentum, origin, b_field, charge - Helix(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q); - - /// ReconstructParticle, b field - Helix(const edm4eic::ReconstructedParticle& p, const double b_field); - - ~Helix() = default; - - double dipAngle() const; - double curvature() const; /// 1/R in xy-plane - double phase() const; /// aziumth in xy-plane measured from ring center - double xcenter() const; /// x-center of circle in xy-plane - double ycenter() const; /// y-center of circle in xy-plane - int h() const; /// -sign(q*B); - - const edm4hep::Vector3f& origin() const; /// starting point - - void setParameters(double c, double dip, double phase, const edm4hep::Vector3f& o, int h); - - void setParameters(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q); - - /// edm4eic::TrackParameters, b field - void setParameters(const edm4eic::TrackParameters& trk, const double b_field); - - /// coordinates of helix at point s - double x(double s) const; - double y(double s) const; - double z(double s) const; - - edm4hep::Vector3f at(double s) const; - - /// pointing vector of helix at point s - double cx(double s) const; - double cy(double s) const; - double cz(double s = 0) const; - - edm4hep::Vector3f cat(double s) const; - - /// returns period length of helix - double period() const; - - /// path length at given r (cylindrical r) - std::pair pathLength(double r) const; - - /// path length at given r (cylindrical r, cylinder axis at x,y) - std::pair pathLength(double r, double x, double y); - - /// path length at distance of closest approach to a given point - double pathLength(const edm4hep::Vector3f& p, bool scanPeriods = true) const; - - /// path length at intersection with plane - double pathLength(const edm4hep::Vector3f& r, - const edm4hep::Vector3f& n) const; - - /// path length at distance of closest approach in the xy-plane to a given point - double pathLength(double x, double y) const; - - /// path lengths at dca between two helices - std::pair pathLengths(const Helix&, - double minStepSize = 10*edm4eic::unit::um, - double minRange = 10*edm4eic::unit::cm) const; - - /// minimal distance between point and helix - double distance(const edm4hep::Vector3f& p, bool scanPeriods = true) const; - - /// checks for valid parametrization - bool valid(double world = 1.e+5) const {return !bad(world);} - int bad(double world = 1.e+5) const; - - /// move the origin along the helix to s which becomes then s=0 - virtual void moveOrigin(double s); - - static const double NoSolution; - - - void setCurvature(double); /// performs also various checks - void setPhase(double); - void setDipAngle(double); - - /// value of S where distance in x-y plane is minimal - double fudgePathLength(const edm4hep::Vector3f&) const; - - // Requires: signed Magnetic Field - edm4hep::Vector3f momentum(double) const; // returns the momentum at origin - edm4hep::Vector3f momentumAt(double, double) const; // returns momemtum at S - int charge(double) const; // returns charge of particle - // 2d DCA to x,y point signed relative to curvature - double curvatureSignedDistance(double x, double y) ; - // 2d DCA to x,y point signed relative to rotation - double geometricSignedDistance(double x, double y) ; - // 3d DCA to 3d point signed relative to curvature - double curvatureSignedDistance(const edm4hep::Vector3f&) ; - // 3d DCA to 3d point signed relative to rotation - double geometricSignedDistance(const edm4hep::Vector3f&) ; - - // - void Print() const; + /// curvature, dip angle, phase, origin, h + Helix(const double c, const double dip, const double phase, const edm4hep::Vector3f& o, + const int h = -1); + + /// momentum, origin, b_field, charge + Helix(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q); + + /// ReconstructParticle, b field + Helix(const edm4eic::ReconstructedParticle& p, const double b_field); + + ~Helix() = default; + + double dipAngle() const; + double curvature() const; /// 1/R in xy-plane + double phase() const; /// aziumth in xy-plane measured from ring center + double xcenter() const; /// x-center of circle in xy-plane + double ycenter() const; /// y-center of circle in xy-plane + int h() const; /// -sign(q*B); + + const edm4hep::Vector3f& origin() const; /// starting point + + void setParameters(double c, double dip, double phase, const edm4hep::Vector3f& o, int h); + + void setParameters(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, + const int q); + + /// edm4eic::TrackParameters, b field + void setParameters(const edm4eic::TrackParameters& trk, const double b_field); + + /// coordinates of helix at point s + double x(double s) const; + double y(double s) const; + double z(double s) const; + + edm4hep::Vector3f at(double s) const; + + /// pointing vector of helix at point s + double cx(double s) const; + double cy(double s) const; + double cz(double s = 0) const; + + edm4hep::Vector3f cat(double s) const; + + /// returns period length of helix + double period() const; + + /// path length at given r (cylindrical r) + std::pair pathLength(double r) const; + + /// path length at given r (cylindrical r, cylinder axis at x,y) + std::pair pathLength(double r, double x, double y); + + /// path length at distance of closest approach to a given point + double pathLength(const edm4hep::Vector3f& p, bool scanPeriods = true) const; + + /// path length at intersection with plane + double pathLength(const edm4hep::Vector3f& r, const edm4hep::Vector3f& n) const; + + /// path length at distance of closest approach in the xy-plane to a given point + double pathLength(double x, double y) const; + + /// path lengths at dca between two helices + std::pair pathLengths(const Helix&, double minStepSize = 10 * edm4eic::unit::um, + double minRange = 10 * edm4eic::unit::cm) const; + + /// minimal distance between point and helix + double distance(const edm4hep::Vector3f& p, bool scanPeriods = true) const; + + /// checks for valid parametrization + bool valid(double world = 1.e+5) const { return !bad(world); } + int bad(double world = 1.e+5) const; + + /// move the origin along the helix to s which becomes then s=0 + virtual void moveOrigin(double s); + + static const double NoSolution; + + void setCurvature(double); /// performs also various checks + void setPhase(double); + void setDipAngle(double); + + /// value of S where distance in x-y plane is minimal + double fudgePathLength(const edm4hep::Vector3f&) const; + + // Requires: signed Magnetic Field + edm4hep::Vector3f momentum(double) const; // returns the momentum at origin + edm4hep::Vector3f momentumAt(double, double) const; // returns momemtum at S + int charge(double) const; // returns charge of particle + // 2d DCA to x,y point signed relative to curvature + double curvatureSignedDistance(double x, double y); + // 2d DCA to x,y point signed relative to rotation + double geometricSignedDistance(double x, double y); + // 3d DCA to 3d point signed relative to curvature + double curvatureSignedDistance(const edm4hep::Vector3f&); + // 3d DCA to 3d point signed relative to rotation + double geometricSignedDistance(const edm4hep::Vector3f&); + + // + void Print() const; }; // end class Helix @@ -147,92 +146,83 @@ class Helix { // // Inline functions // -inline int Helix::h() const {return mH;} +inline int Helix::h() const { return mH; } -inline double Helix::dipAngle() const {return mDipAngle;} +inline double Helix::dipAngle() const { return mDipAngle; } -inline double Helix::curvature() const {return mCurvature;} +inline double Helix::curvature() const { return mCurvature; } -inline double Helix::phase() const {return mPhase;} +inline double Helix::phase() const { return mPhase; } -inline double Helix::x(double s) const -{ - if (mSingularity) - return mOrigin.x - s*mCosDipAngle*mSinPhase; - else - return mOrigin.x + (cos(mPhase + s*mH*mCurvature*mCosDipAngle)-mCosPhase)/mCurvature; -} - -inline double Helix::y(double s) const -{ - if (mSingularity) - return mOrigin.y + s*mCosDipAngle*mCosPhase; - else - return mOrigin.y + (sin(mPhase + s*mH*mCurvature*mCosDipAngle)-mSinPhase)/mCurvature; +inline double Helix::x(double s) const { + if (mSingularity) + return mOrigin.x - s * mCosDipAngle * mSinPhase; + else + return mOrigin.x + (cos(mPhase + s * mH * mCurvature * mCosDipAngle) - mCosPhase) / mCurvature; } -inline double Helix::z(double s) const -{ - return mOrigin.z + s*mSinDipAngle; +inline double Helix::y(double s) const { + if (mSingularity) + return mOrigin.y + s * mCosDipAngle * mCosPhase; + else + return mOrigin.y + (sin(mPhase + s * mH * mCurvature * mCosDipAngle) - mSinPhase) / mCurvature; } -inline double Helix::cx(double s) const -{ - if (mSingularity) - return -mCosDipAngle*mSinPhase; - else - return -sin(mPhase + s*mH*mCurvature*mCosDipAngle)*mH*mCosDipAngle; +inline double Helix::z(double s) const { return mOrigin.z + s * mSinDipAngle; } + +inline double Helix::cx(double s) const { + if (mSingularity) + return -mCosDipAngle * mSinPhase; + else + return -sin(mPhase + s * mH * mCurvature * mCosDipAngle) * mH * mCosDipAngle; } -inline double Helix::cy(double s) const -{ - if (mSingularity) - return mCosDipAngle*mCosPhase; - else - return cos(mPhase + s*mH*mCurvature*mCosDipAngle)*mH*mCosDipAngle; +inline double Helix::cy(double s) const { + if (mSingularity) + return mCosDipAngle * mCosPhase; + else + return cos(mPhase + s * mH * mCurvature * mCosDipAngle) * mH * mCosDipAngle; } -inline double Helix::cz(double /* s */) const -{ - return mSinDipAngle; -} +inline double Helix::cz(double /* s */) const { return mSinDipAngle; } -inline const edm4hep::Vector3f& Helix::origin() const {return mOrigin;} +inline const edm4hep::Vector3f& Helix::origin() const { return mOrigin; } -inline edm4hep::Vector3f Helix::at(double s) const -{ - return edm4hep::Vector3f(x(s), y(s), z(s)); -} +inline edm4hep::Vector3f Helix::at(double s) const { return edm4hep::Vector3f(x(s), y(s), z(s)); } -inline edm4hep::Vector3f Helix::cat(double s) const -{ - return edm4hep::Vector3f(cx(s), cy(s), cz(s)); +inline edm4hep::Vector3f Helix::cat(double s) const { + return edm4hep::Vector3f(cx(s), cy(s), cz(s)); } -inline double Helix::pathLength(double X, double Y) const -{ - return fudgePathLength(edm4hep::Vector3f(X, Y, 0)); +inline double Helix::pathLength(double X, double Y) const { + return fudgePathLength(edm4hep::Vector3f(X, Y, 0)); } -inline int Helix::bad(double WorldSize) const -{ +inline int Helix::bad(double WorldSize) const { -// int ierr; - if (!::finite(mDipAngle )) return 11; - if (!::finite(mCurvature )) return 12; + // int ierr; + if (!::finite(mDipAngle)) + return 11; + if (!::finite(mCurvature)) + return 12; -// ierr = mOrigin.bad(WorldSize); -// if (ierr) return 3+ierr*100; + // ierr = mOrigin.bad(WorldSize); + // if (ierr) return 3+ierr*100; - if (::fabs(mDipAngle) >1.58) return 21; - double qwe = ::fabs(::fabs(mDipAngle)-M_PI/2); - if (qwe < 1./WorldSize ) return 31; + if (::fabs(mDipAngle) > 1.58) + return 21; + double qwe = ::fabs(::fabs(mDipAngle) - M_PI / 2); + if (qwe < 1. / WorldSize) + return 31; - if (::fabs(mCurvature) > WorldSize) return 22; - if (mCurvature < 0 ) return 32; + if (::fabs(mCurvature) > WorldSize) + return 22; + if (mCurvature < 0) + return 32; - if (abs(mH) != 1 ) return 24; + if (abs(mH) != 1) + return 24; - return 0; + return 0; } -} +} // namespace eicrecon diff --git a/src/algorithms/reco/SecondaryVerticesHelix.cc b/src/algorithms/reco/SecondaryVerticesHelix.cc index e1fc0d2ec7..40c65e4c8e 100644 --- a/src/algorithms/reco/SecondaryVerticesHelix.cc +++ b/src/algorithms/reco/SecondaryVerticesHelix.cc @@ -36,101 +36,115 @@ void SecondaryVerticesHelix::init() {} */ void SecondaryVerticesHelix::process(const SecondaryVerticesHelix::Input& input, const SecondaryVerticesHelix::Output& output) const { - const auto [rcvtx, rcparts] = input; + const auto [rcvtx, rcparts] = input; auto [out_secondary_vertices] = output; auto& particleSvc = algorithms::ParticleSvc::instance(); - if((*rcvtx).size()==0) { - info(" No primary vertex in this event! Skip secondary vertex finder!"); + if ((*rcvtx).size() == 0) { + info(" No primary vertex in this event! Skip secondary vertex finder!"); return; } const auto pVtxPos4f = (*rcvtx)[0].getPosition(); // convert to cm - edm4hep::Vector3f pVtxPos(pVtxPos4f.x*edm4eic::unit::mm/edm4eic::unit::cm, - pVtxPos4f.y*edm4eic::unit::mm/edm4eic::unit::cm, - pVtxPos4f.z*edm4eic::unit::mm/edm4eic::unit::cm); - info("\t Primary vertex = ({},{},{})cm \t b field = {} tesla", pVtxPos.x, pVtxPos.y, pVtxPos.z, m_cfg.b_field/dd4hep::tesla); + edm4hep::Vector3f pVtxPos(pVtxPos4f.x * edm4eic::unit::mm / edm4eic::unit::cm, + pVtxPos4f.y * edm4eic::unit::mm / edm4eic::unit::cm, + pVtxPos4f.z * edm4eic::unit::mm / edm4eic::unit::cm); + info("\t Primary vertex = ({},{},{})cm \t b field = {} tesla", pVtxPos.x, pVtxPos.y, pVtxPos.z, + m_cfg.b_field / dd4hep::tesla); - std::vector hVec; hVec.clear(); - std::vector indexVec; indexVec.clear(); + std::vector hVec; + hVec.clear(); + std::vector indexVec; + indexVec.clear(); for (unsigned int i = 0; const auto& p : *rcparts) { - if(p.getCharge()==0) continue; + if (p.getCharge() == 0) + continue; Helix h(p, m_cfg.b_field); double dca = h.distance(pVtxPos) * edm4eic::unit::cm; - if( dca < m_cfg.minDca ) continue; - + if (dca < m_cfg.minDca) + continue; + hVec.push_back(h); indexVec.push_back(i); ++i; } - - if( hVec.size() != indexVec.size() ) return; - + + if (hVec.size() != indexVec.size()) + return; + debug("\t Vector size {}, {}", hVec.size(), indexVec.size()); for (unsigned int i1 = 0; i1 < hVec.size(); ++i1) { for (unsigned int i2 = i1 + 1; i2 < hVec.size(); ++i2) { - const auto& p1 = (*rcparts)[indexVec[i1]]; - const auto& p2 = (*rcparts)[indexVec[i2]]; - - if ( !(m_cfg.unlikesign && p1.getCharge() + p2.getCharge() == 0) ) continue; - - const auto& h1 = hVec[i1]; - const auto& h2 = hVec[i2]; - - // Helix function uses cm unit - double dca1 = h1.distance(pVtxPos) * edm4eic::unit::cm; - double dca2 = h2.distance(pVtxPos) * edm4eic::unit::cm; - if( dca1 < m_cfg.minDca || dca2 < m_cfg.minDca ) continue; - - std::pair const ss = h1.pathLengths(h2); - edm4hep::Vector3f h1AtDcaTo2 = h1.at(ss.first); - edm4hep::Vector3f h2AtDcaTo1 = h2.at(ss.second); - - double dca12 = edm4hep::utils::magnitude(h1AtDcaTo2 - h2AtDcaTo1) * edm4eic::unit::cm; - if(std::isnan(dca12)) continue; - if( dca12 > m_cfg.maxDca12 ) continue; - edm4hep::Vector3f pairPos = 0.5*(h1AtDcaTo2 + h2AtDcaTo1); - - edm4hep::Vector3f h1MomAtDca = h1.momentumAt(ss.first, m_cfg.b_field); - edm4hep::Vector3f h2MomAtDca = h2.momentumAt(ss.second, m_cfg.b_field); - edm4hep::Vector3f pairMom = h1MomAtDca + h2MomAtDca; - - double e1 = std::hypot(edm4hep::utils::magnitude(h1MomAtDca), particleSvc.particle(p1.getPDG()).mass); - double e2 = std::hypot(edm4hep::utils::magnitude(h2MomAtDca), particleSvc.particle(p2.getPDG()).mass); - double pairE = e1+e2; - double pairP = edm4hep::utils::magnitude(pairMom); - - double m_inv2 = pairE*pairE - pairP*pairP; - double m_inv = (m_inv2>0) ? sqrt(m_inv2) : 0.; - double angle = edm4hep::utils::angleBetween(pairMom, pairPos - pVtxPos); - if(cos(angle) < m_cfg.minCostheta ) continue; - - double beta = edm4hep::utils::magnitude(pairMom)/pairE; - double time = edm4hep::utils::magnitude(pairPos - pVtxPos)/(beta*dd4hep::c_light); - edm4hep::Vector3f dL = pairPos - pVtxPos; // in cm - edm4hep::Vector3f decayL(dL.x*edm4eic::unit::cm, dL.y*edm4eic::unit::cm, dL.z*edm4eic::unit::cm); - double dca2pv = edm4hep::utils::magnitude(decayL) * sin(angle); - if(dca2pv > m_cfg.maxDca) continue; - - auto v0 = out_secondary_vertices->create(); - v0.setType(2); // 2 for secondary - v0.setPosition({(float)(pairPos.x * edm4eic::unit::cm / edm4eic::unit::mm), - (float)(pairPos.y * edm4eic::unit::cm / edm4eic::unit::mm), - (float)(pairPos.z * edm4eic::unit::cm / edm4eic::unit::mm), - (float)time}); - v0.addToAssociatedParticles(p1); - v0.addToAssociatedParticles(p2); - - info("One secondary vertex found at (x,y,z) = ({}, {}, {}) mm.", - pairPos.x * edm4eic::unit::cm / edm4eic::unit::mm, - pairPos.y * edm4eic::unit::cm / edm4eic::unit::mm, - pairPos.x * edm4eic::unit::cm / edm4eic::unit::mm); - + const auto& p1 = (*rcparts)[indexVec[i1]]; + const auto& p2 = (*rcparts)[indexVec[i2]]; + + if (!(m_cfg.unlikesign && p1.getCharge() + p2.getCharge() == 0)) + continue; + + const auto& h1 = hVec[i1]; + const auto& h2 = hVec[i2]; + + // Helix function uses cm unit + double dca1 = h1.distance(pVtxPos) * edm4eic::unit::cm; + double dca2 = h2.distance(pVtxPos) * edm4eic::unit::cm; + if (dca1 < m_cfg.minDca || dca2 < m_cfg.minDca) + continue; + + std::pair const ss = h1.pathLengths(h2); + edm4hep::Vector3f h1AtDcaTo2 = h1.at(ss.first); + edm4hep::Vector3f h2AtDcaTo1 = h2.at(ss.second); + + double dca12 = edm4hep::utils::magnitude(h1AtDcaTo2 - h2AtDcaTo1) * edm4eic::unit::cm; + if (std::isnan(dca12)) + continue; + if (dca12 > m_cfg.maxDca12) + continue; + edm4hep::Vector3f pairPos = 0.5 * (h1AtDcaTo2 + h2AtDcaTo1); + + edm4hep::Vector3f h1MomAtDca = h1.momentumAt(ss.first, m_cfg.b_field); + edm4hep::Vector3f h2MomAtDca = h2.momentumAt(ss.second, m_cfg.b_field); + edm4hep::Vector3f pairMom = h1MomAtDca + h2MomAtDca; + + double e1 = + std::hypot(edm4hep::utils::magnitude(h1MomAtDca), particleSvc.particle(p1.getPDG()).mass); + double e2 = + std::hypot(edm4hep::utils::magnitude(h2MomAtDca), particleSvc.particle(p2.getPDG()).mass); + double pairE = e1 + e2; + double pairP = edm4hep::utils::magnitude(pairMom); + + double m_inv2 = pairE * pairE - pairP * pairP; + double m_inv = (m_inv2 > 0) ? sqrt(m_inv2) : 0.; + double angle = edm4hep::utils::angleBetween(pairMom, pairPos - pVtxPos); + if (cos(angle) < m_cfg.minCostheta) + continue; + + double beta = edm4hep::utils::magnitude(pairMom) / pairE; + double time = edm4hep::utils::magnitude(pairPos - pVtxPos) / (beta * dd4hep::c_light); + edm4hep::Vector3f dL = pairPos - pVtxPos; // in cm + edm4hep::Vector3f decayL(dL.x * edm4eic::unit::cm, dL.y * edm4eic::unit::cm, + dL.z * edm4eic::unit::cm); + double dca2pv = edm4hep::utils::magnitude(decayL) * sin(angle); + if (dca2pv > m_cfg.maxDca) + continue; + + auto v0 = out_secondary_vertices->create(); + v0.setType(2); // 2 for secondary + v0.setPosition({(float)(pairPos.x * edm4eic::unit::cm / edm4eic::unit::mm), + (float)(pairPos.y * edm4eic::unit::cm / edm4eic::unit::mm), + (float)(pairPos.z * edm4eic::unit::cm / edm4eic::unit::mm), (float)time}); + v0.addToAssociatedParticles(p1); + v0.addToAssociatedParticles(p2); + + info("One secondary vertex found at (x,y,z) = ({}, {}, {}) mm.", + pairPos.x * edm4eic::unit::cm / edm4eic::unit::mm, + pairPos.y * edm4eic::unit::cm / edm4eic::unit::mm, + pairPos.x * edm4eic::unit::cm / edm4eic::unit::mm); + } // end i2 } // end i1 - + } // end process } // namespace eicrecon diff --git a/src/algorithms/reco/SecondaryVerticesHelix.h b/src/algorithms/reco/SecondaryVerticesHelix.h index dc3ece37a2..cf065d2ec6 100644 --- a/src/algorithms/reco/SecondaryVerticesHelix.h +++ b/src/algorithms/reco/SecondaryVerticesHelix.h @@ -19,14 +19,15 @@ using SecondaryVerticesHelixAlgorithm = algorithms::Algorithm< algorithms::Output>; class SecondaryVerticesHelix : public SecondaryVerticesHelixAlgorithm, - public WithPodConfig { + public WithPodConfig { public: SecondaryVerticesHelix(std::string_view name) - : SecondaryVerticesHelixAlgorithm{name, - {"inputVertices", "inputParticles"}, - {"outputSecondaryVertices"}, - "Reconstruct secondary vertices in SecondaryVertices collection"} {} + : SecondaryVerticesHelixAlgorithm{ + name, + {"inputVertices", "inputParticles"}, + {"outputSecondaryVertices"}, + "Reconstruct secondary vertices in SecondaryVertices collection"} {} void init() final; void process(const Input&, const Output&) const final; diff --git a/src/algorithms/reco/SecondaryVerticesHelixConfig.h b/src/algorithms/reco/SecondaryVerticesHelixConfig.h index 3321e9fab5..cef2d8e300 100644 --- a/src/algorithms/reco/SecondaryVerticesHelixConfig.h +++ b/src/algorithms/reco/SecondaryVerticesHelixConfig.h @@ -11,11 +11,11 @@ namespace eicrecon { struct SecondaryVerticesHelixConfig { - float b_field = -1.7*dd4hep::tesla; - bool unlikesign = true; - float minDca = 0.03*edm4eic::unit::mm; // mm, daughter to pVtx - float maxDca12 = 1.*edm4eic::unit::mm; // mm, dca between daughter 1 and 2 - float maxDca = 1.*edm4eic::unit::mm; // mm, dca of V0 to pVtx + float b_field = -1.7 * dd4hep::tesla; + bool unlikesign = true; + float minDca = 0.03 * edm4eic::unit::mm; // mm, daughter to pVtx + float maxDca12 = 1. * edm4eic::unit::mm; // mm, dca between daughter 1 and 2 + float maxDca = 1. * edm4eic::unit::mm; // mm, dca of V0 to pVtx float minCostheta = 0.8; // costheta, theta: angle of V0 decay direction and momentum }; diff --git a/src/factories/reco/SecondaryVerticesHelix_factory.h b/src/factories/reco/SecondaryVerticesHelix_factory.h index 17fa2df7b5..b52a459272 100644 --- a/src/factories/reco/SecondaryVerticesHelix_factory.h +++ b/src/factories/reco/SecondaryVerticesHelix_factory.h @@ -43,7 +43,8 @@ class SecondaryVerticesHelix_factory } void Process(int32_t /* run_number */, uint64_t /* event_number */) { - m_algo->process({m_rc_vertices_input(), m_rc_parts_input()}, {m_secondary_vertices_output().get()}); + m_algo->process({m_rc_vertices_input(), m_rc_parts_input()}, + {m_secondary_vertices_output().get()}); } }; diff --git a/src/global/reco/reco.cc b/src/global/reco/reco.cc index a934afec99..51cd0af995 100644 --- a/src/global/reco/reco.cc +++ b/src/global/reco/reco.cc @@ -276,7 +276,7 @@ void InitPlugin(JApplication* app) { "PrimaryVertices", {"CentralTrackVertices"}, {"PrimaryVertices"}, {}, app)); app->Add(new JOmniFactoryGeneratorT( - "SecondaryVerticesHelix", {"PrimaryVertices","ReconstructedParticles"}, {"SecondaryVerticesHelix"}, {}, app)); - + "SecondaryVerticesHelix", {"PrimaryVertices", "ReconstructedParticles"}, + {"SecondaryVerticesHelix"}, {}, app)); } } // extern "C" From fe5f71e46ff244fe0f864c845f6ab5a19ac22839 Mon Sep 17 00:00:00 2001 From: Xin Dong Date: Mon, 20 Oct 2025 17:24:33 -0400 Subject: [PATCH 08/17] fixed typos in Helix function (in commenting area) --- src/algorithms/reco/Helix.cc | 125 +++++++++++++++++++++++++++++++++++ src/algorithms/reco/Helix.h | 18 +++++ 2 files changed, 143 insertions(+) diff --git a/src/algorithms/reco/Helix.cc b/src/algorithms/reco/Helix.cc index 6664fb3f8e..9b1338cb93 100644 --- a/src/algorithms/reco/Helix.cc +++ b/src/algorithms/reco/Helix.cc @@ -514,6 +514,7 @@ std::pair Helix::pathLengths(const Helix& h, double minStepSize, } } +<<<<<<< HEAD void Helix::moveOrigin(double s) { if (mSingularity) mOrigin = at(s); @@ -523,6 +524,130 @@ void Helix::moveOrigin(double s) { mOrigin = newOrigin; setPhase(newPhase); } +======= +std::pair +Helix::pathLengths(const Helix& h, double minStepSize, double minRange) const +{ + // + // Cannot handle case where one is a helix + // and the other one is a straight line. + // + if (mSingularity != h.mSingularity) + return std::pair(NoSolution, NoSolution); + + double s1, s2; + + if (mSingularity) { + // + // Analytic solution + // + edm4hep::Vector3f dv = h.mOrigin - mOrigin; + edm4hep::Vector3f a(-mCosDipAngle*mSinPhase, + mCosDipAngle*mCosPhase, + mSinDipAngle); + edm4hep::Vector3f b(-h.mCosDipAngle*h.mSinPhase, + h.mCosDipAngle*h.mCosPhase, + h.mSinDipAngle); + double ab = a*b; + double g = dv*a; + double k = dv*b; + s2 = (k-ab*g)/(ab*ab-1.); + s1 = g+s2*ab; + return std::pair(s1, s2); + } + else { + // + // First step: get dca in the xy-plane as start value + // + double dx = h.xcenter() - xcenter(); + double dy = h.ycenter() - ycenter(); + double dd = ::sqrt(dx*dx + dy*dy); + double r1 = 1/curvature(); + double r2 = 1/h.curvature(); + + double cosAlpha = (r1*r1 + dd*dd - r2*r2)/(2*r1*dd); + + double s; + double x, y; + if (fabs(cosAlpha) < 1) { // two solutions + double sinAlpha = sin(acos(cosAlpha)); + x = xcenter() + r1*(cosAlpha*dx - sinAlpha*dy)/dd; + y = ycenter() + r1*(sinAlpha*dx + cosAlpha*dy)/dd; + s = pathLength(x, y); + x = xcenter() + r1*(cosAlpha*dx + sinAlpha*dy)/dd; + y = ycenter() + r1*(cosAlpha*dy - sinAlpha*dx)/dd; + double a = pathLength(x, y); + if (h.distance(at(a)) < h.distance(at(s))) s = a; + } + else { // no intersection (or exactly one) + int rsign = ((r2-r1) > dd ? -1 : 1); // set -1 when *this* helix is + // completely contained in the other + x = xcenter() + rsign*r1*dx/dd; + y = ycenter() + rsign*r1*dy/dd; + s = pathLength(x, y); + } + + // + // Second step: scan in decreasing intervals around seed 's' + // minRange and minStepSize are passed as arguments to the method. + // They have default values defined in the header file. + // + double dmin = h.distance(at(s)); + double range = std::max(2*dmin, minRange); + double ds = range/10; + double slast=-999999, ss, d; + s1 = s - range/2.; + s2 = s + range/2.; + + while (ds > minStepSize) { + for (ss=s1; ss(s, h.pathLength(at(s))); + } +} + + +void Helix::moveOrigin(double s) +{ + if (mSingularity) + mOrigin = at(s); + else { + edm4hep::Vector3f newOrigin = at(s); + double newPhase = atan2(newOrigin.y - ycenter(), + newOrigin.x - xcenter()); + mOrigin = newOrigin; + setPhase(newPhase); + } +>>>>>>> 5c2ad0d0 (fixed typos in Helix function (in commenting area)) } /* int operator== (const Helix& a, const Helix& b) diff --git a/src/algorithms/reco/Helix.h b/src/algorithms/reco/Helix.h index d68354b7b1..dab146698b 100644 --- a/src/algorithms/reco/Helix.h +++ b/src/algorithms/reco/Helix.h @@ -77,6 +77,7 @@ class Helix { edm4hep::Vector3f cat(double s) const; +<<<<<<< HEAD /// returns period length of helix double period() const; @@ -133,6 +134,23 @@ class Helix { // void Print() const; +======= + // Requires: signed Magnetic Field + edm4hep::Vector3f momentum(double) const; // returns the momentum at origin + edm4hep::Vector3f momentumAt(double, double) const; // returns momentum at S + int charge(double) const; // returns charge of particle + // 2d DCA to x,y point signed relative to curvature + double curvatureSignedDistance(double x, double y) ; + // 2d DCA to x,y point signed relative to rotation + double geometricSignedDistance(double x, double y) ; + // 3d DCA to 3d point signed relative to curvature + double curvatureSignedDistance(const edm4hep::Vector3f&) ; + // 3d DCA to 3d point signed relative to rotation + double geometricSignedDistance(const edm4hep::Vector3f&) ; + + // + void Print() const; +>>>>>>> 5c2ad0d0 (fixed typos in Helix function (in commenting area)) }; // end class Helix From 18b043483f92f9419a257b75e8481abca6a0824f Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 20 Oct 2025 21:30:58 +0000 Subject: [PATCH 09/17] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/algorithms/reco/Helix.cc | 220 +++++++++++++++++------------------ src/algorithms/reco/Helix.h | 30 ++--- 2 files changed, 120 insertions(+), 130 deletions(-) diff --git a/src/algorithms/reco/Helix.cc b/src/algorithms/reco/Helix.cc index 9b1338cb93..46d07c74b1 100644 --- a/src/algorithms/reco/Helix.cc +++ b/src/algorithms/reco/Helix.cc @@ -525,128 +525,118 @@ void Helix::moveOrigin(double s) { setPhase(newPhase); } ======= -std::pair -Helix::pathLengths(const Helix& h, double minStepSize, double minRange) const -{ +std::pair Helix::pathLengths(const Helix& h, double minStepSize, + double minRange) const { + // + // Cannot handle case where one is a helix + // and the other one is a straight line. + // + if (mSingularity != h.mSingularity) + return std::pair(NoSolution, NoSolution); + + double s1, s2; + + if (mSingularity) { + // + // Analytic solution + // + edm4hep::Vector3f dv = h.mOrigin - mOrigin; + edm4hep::Vector3f a(-mCosDipAngle * mSinPhase, mCosDipAngle * mCosPhase, mSinDipAngle); + edm4hep::Vector3f b(-h.mCosDipAngle * h.mSinPhase, h.mCosDipAngle * h.mCosPhase, + h.mSinDipAngle); + double ab = a * b; + double g = dv * a; + double k = dv * b; + s2 = (k - ab * g) / (ab * ab - 1.); + s1 = g + s2 * ab; + return std::pair(s1, s2); + } else { // - // Cannot handle case where one is a helix - // and the other one is a straight line. + // First step: get dca in the xy-plane as start value // - if (mSingularity != h.mSingularity) - return std::pair(NoSolution, NoSolution); - - double s1, s2; - - if (mSingularity) { - // - // Analytic solution - // - edm4hep::Vector3f dv = h.mOrigin - mOrigin; - edm4hep::Vector3f a(-mCosDipAngle*mSinPhase, - mCosDipAngle*mCosPhase, - mSinDipAngle); - edm4hep::Vector3f b(-h.mCosDipAngle*h.mSinPhase, - h.mCosDipAngle*h.mCosPhase, - h.mSinDipAngle); - double ab = a*b; - double g = dv*a; - double k = dv*b; - s2 = (k-ab*g)/(ab*ab-1.); - s1 = g+s2*ab; - return std::pair(s1, s2); + double dx = h.xcenter() - xcenter(); + double dy = h.ycenter() - ycenter(); + double dd = ::sqrt(dx * dx + dy * dy); + double r1 = 1 / curvature(); + double r2 = 1 / h.curvature(); + + double cosAlpha = (r1 * r1 + dd * dd - r2 * r2) / (2 * r1 * dd); + + double s; + double x, y; + if (fabs(cosAlpha) < 1) { // two solutions + double sinAlpha = sin(acos(cosAlpha)); + x = xcenter() + r1 * (cosAlpha * dx - sinAlpha * dy) / dd; + y = ycenter() + r1 * (sinAlpha * dx + cosAlpha * dy) / dd; + s = pathLength(x, y); + x = xcenter() + r1 * (cosAlpha * dx + sinAlpha * dy) / dd; + y = ycenter() + r1 * (cosAlpha * dy - sinAlpha * dx) / dd; + double a = pathLength(x, y); + if (h.distance(at(a)) < h.distance(at(s))) + s = a; + } else { // no intersection (or exactly one) + int rsign = ((r2 - r1) > dd ? -1 : 1); // set -1 when *this* helix is + // completely contained in the other + x = xcenter() + rsign * r1 * dx / dd; + y = ycenter() + rsign * r1 * dy / dd; + s = pathLength(x, y); } - else { - // - // First step: get dca in the xy-plane as start value - // - double dx = h.xcenter() - xcenter(); - double dy = h.ycenter() - ycenter(); - double dd = ::sqrt(dx*dx + dy*dy); - double r1 = 1/curvature(); - double r2 = 1/h.curvature(); - - double cosAlpha = (r1*r1 + dd*dd - r2*r2)/(2*r1*dd); - - double s; - double x, y; - if (fabs(cosAlpha) < 1) { // two solutions - double sinAlpha = sin(acos(cosAlpha)); - x = xcenter() + r1*(cosAlpha*dx - sinAlpha*dy)/dd; - y = ycenter() + r1*(sinAlpha*dx + cosAlpha*dy)/dd; - s = pathLength(x, y); - x = xcenter() + r1*(cosAlpha*dx + sinAlpha*dy)/dd; - y = ycenter() + r1*(cosAlpha*dy - sinAlpha*dx)/dd; - double a = pathLength(x, y); - if (h.distance(at(a)) < h.distance(at(s))) s = a; - } - else { // no intersection (or exactly one) - int rsign = ((r2-r1) > dd ? -1 : 1); // set -1 when *this* helix is - // completely contained in the other - x = xcenter() + rsign*r1*dx/dd; - y = ycenter() + rsign*r1*dy/dd; - s = pathLength(x, y); - } - - // - // Second step: scan in decreasing intervals around seed 's' - // minRange and minStepSize are passed as arguments to the method. - // They have default values defined in the header file. - // - double dmin = h.distance(at(s)); - double range = std::max(2*dmin, minRange); - double ds = range/10; - double slast=-999999, ss, d; - s1 = s - range/2.; - s2 = s + range/2.; - - while (ds > minStepSize) { - for (ss=s1; ss minStepSize) { + for (ss = s1; ss < s2 + ds; ss += ds) { + d = h.distance(at(ss)); + if (d < dmin) { + dmin = d; + s = ss; } - return std::pair(s, h.pathLength(at(s))); + slast = ss; + } + // + // In the rare cases where the minimum is at the + // the border of the current range we shift the range + // and start all over, i.e we do not decrease 'ds'. + // Else we decrease the search interval around the + // current minimum and redo the scan in smaller steps. + // + if (s == s1) { + d = 0.8 * (s2 - s1); + s1 -= d; + s2 -= d; + } else if (s == slast) { + d = 0.8 * (s2 - s1); + s1 += d; + s2 += d; + } else { + s1 = s - ds; + s2 = s + ds; + ds /= 10; + } } + return std::pair(s, h.pathLength(at(s))); + } } - -void Helix::moveOrigin(double s) -{ - if (mSingularity) - mOrigin = at(s); - else { - edm4hep::Vector3f newOrigin = at(s); - double newPhase = atan2(newOrigin.y - ycenter(), - newOrigin.x - xcenter()); - mOrigin = newOrigin; - setPhase(newPhase); - } +void Helix::moveOrigin(double s) { + if (mSingularity) + mOrigin = at(s); + else { + edm4hep::Vector3f newOrigin = at(s); + double newPhase = atan2(newOrigin.y - ycenter(), newOrigin.x - xcenter()); + mOrigin = newOrigin; + setPhase(newPhase); + } >>>>>>> 5c2ad0d0 (fixed typos in Helix function (in commenting area)) } /* diff --git a/src/algorithms/reco/Helix.h b/src/algorithms/reco/Helix.h index dab146698b..2f30c1d6ce 100644 --- a/src/algorithms/reco/Helix.h +++ b/src/algorithms/reco/Helix.h @@ -135,21 +135,21 @@ class Helix { // void Print() const; ======= - // Requires: signed Magnetic Field - edm4hep::Vector3f momentum(double) const; // returns the momentum at origin - edm4hep::Vector3f momentumAt(double, double) const; // returns momentum at S - int charge(double) const; // returns charge of particle - // 2d DCA to x,y point signed relative to curvature - double curvatureSignedDistance(double x, double y) ; - // 2d DCA to x,y point signed relative to rotation - double geometricSignedDistance(double x, double y) ; - // 3d DCA to 3d point signed relative to curvature - double curvatureSignedDistance(const edm4hep::Vector3f&) ; - // 3d DCA to 3d point signed relative to rotation - double geometricSignedDistance(const edm4hep::Vector3f&) ; - - // - void Print() const; + // Requires: signed Magnetic Field + edm4hep::Vector3f momentum(double) const; // returns the momentum at origin + edm4hep::Vector3f momentumAt(double, double) const; // returns momentum at S + int charge(double) const; // returns charge of particle + // 2d DCA to x,y point signed relative to curvature + double curvatureSignedDistance(double x, double y); + // 2d DCA to x,y point signed relative to rotation + double geometricSignedDistance(double x, double y); + // 3d DCA to 3d point signed relative to curvature + double curvatureSignedDistance(const edm4hep::Vector3f&); + // 3d DCA to 3d point signed relative to rotation + double geometricSignedDistance(const edm4hep::Vector3f&); + + // + void Print() const; >>>>>>> 5c2ad0d0 (fixed typos in Helix function (in commenting area)) }; // end class Helix From 41bddbccb721a75e59243c44c3645c63cd9377c8 Mon Sep 17 00:00:00 2001 From: Xin Dong Date: Mon, 20 Oct 2025 18:57:21 -0400 Subject: [PATCH 10/17] fix Helix.* file issue in the last commit after rebasing --- src/algorithms/reco/Helix.cc | 1138 ++++++++++++++++------------------ src/algorithms/reco/Helix.h | 356 ++++++----- 2 files changed, 706 insertions(+), 788 deletions(-) diff --git a/src/algorithms/reco/Helix.cc b/src/algorithms/reco/Helix.cc index 46d07c74b1..a675330c1f 100644 --- a/src/algorithms/reco/Helix.cc +++ b/src/algorithms/reco/Helix.cc @@ -25,13 +25,15 @@ namespace eicrecon { const double Helix::NoSolution = 3.e+33; // basic constructor -Helix::Helix(double c, double d, double phase, const edm4hep::Vector3f& o, int h) { - setParameters(c, d, phase, o, h); +Helix::Helix(double c, double d, double phase, const edm4hep::Vector3f& o, int h) +{ + setParameters(c, d, phase, o, h); } // momentum in GeV, position in cm, B in Tesla -Helix::Helix(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q) { - setParameters(p, o, B, q); +Helix::Helix(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q) +{ + setParameters(p, o, B, q); } // construct using ReconstructParticle @@ -48,596 +50,508 @@ Helix::Helix(const edm4eic::ReconstructedParticle& p, const double b_field) { // construct using TrackParameteters void Helix::setParameters(const edm4eic::TrackParameters& trk, const double b_field) { - const auto mom = edm4hep::utils::sphericalToVector(1.0 / std::abs(trk.getQOverP()), - trk.getTheta(), trk.getPhi()); + const auto mom = edm4hep::utils::sphericalToVector(1.0 / std::abs(trk.getQOverP()), trk.getTheta(), trk.getPhi()); const auto charge = std::copysign(1., trk.getQOverP()); - const auto phi = trk.getPhi(); - const auto loc = trk.getLoc(); - edm4hep::Vector3f pos(-1. * loc.a * sin(phi), loc.a * cos(phi), loc.b); // PCA point + const auto phi = trk.getPhi(); + const auto loc = trk.getLoc(); + edm4hep::Vector3f pos( -1. * loc.a * sin(phi), loc.a * cos(phi), loc.b); // PCA point - setParameters(mom * edm4eic::unit::GeV / dd4hep::GeV, pos * edm4eic::unit::mm / edm4eic::unit::cm, - b_field, charge); + setParameters(mom*edm4eic::unit::GeV/dd4hep::GeV, pos*edm4eic::unit::mm/edm4eic::unit::cm, b_field, charge); } -void Helix::setParameters(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, - const int q) { - mH = (q * B <= 0) ? 1 : -1; - if (p.y == 0 && p.x == 0) - setPhase((M_PI / 4) * (1 - 2. * mH)); - else - setPhase(atan2(p.y, p.x) - mH * M_PI / 2); - setDipAngle(atan2(p.z, edm4hep::utils::magnitudeTransverse(p))); - mOrigin = o; - - double curvature_val = - fabs((dd4hep::c_light * dd4hep::nanosecond / dd4hep::meter * q * B / dd4hep::tesla) / - (edm4hep::utils::magnitude(p) / dd4hep::GeV * mCosDipAngle) / dd4hep::meter); - setCurvature(curvature_val); -} +void Helix::setParameters(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q) +{ + mH = (q*B <= 0) ? 1 : -1; + if(p.y == 0 && p.x == 0) + setPhase((M_PI/4)*(1-2.*mH)); + else + setPhase(atan2(p.y,p.x)-mH*M_PI/2); + setDipAngle(atan2(p.z,edm4hep::utils::magnitudeTransverse(p))); + mOrigin = o; -void Helix::setParameters(double c, double dip, double phase, const edm4hep::Vector3f& o, int h) { - // - // The order in which the parameters are set is important - // since setCurvature might have to adjust the others. - // - mH = (h >= 0) ? 1 : -1; // Default is: positive particle - // positive field - mOrigin = o; - setDipAngle(dip); - setPhase(phase); - - // - // Check for singularity and correct for negative curvature. - // May change mH and mPhase. Must therefore be set last. - // - setCurvature(c); - - // - // For the case B=0, h is ill defined. In the following we - // always assume h = +1. Since phase = psi - h * pi/2 - // we have to correct the phase in case h = -1. - // This assumes that the user uses the same h for phase - // as the one he passed to the constructor. - // - if (mSingularity && mH == -1) { - mH = +1; - setPhase(mPhase - M_PI); - } + double curvature_val = fabs((dd4hep::c_light*dd4hep::nanosecond/dd4hep::meter*q*B/dd4hep::tesla)/ + (edm4hep::utils::magnitude(p)/dd4hep::GeV*mCosDipAngle)/dd4hep::meter); + setCurvature(curvature_val); } -void Helix::setCurvature(double val) { - if (val < 0) { - mCurvature = -val; - mH = -mH; - setPhase(mPhase + M_PI); - } else - mCurvature = val; - - if (fabs(mCurvature) <= std::numeric_limits::epsilon()) - mSingularity = true; // straight line - else - mSingularity = false; // curved -} +void Helix::setParameters(double c, double dip, double phase, + const edm4hep::Vector3f& o, int h) +{ + // + // The order in which the parameters are set is important + // since setCurvature might have to adjust the others. + // + mH = (h>=0) ? 1 : -1; // Default is: positive particle + // positive field + mOrigin = o; + setDipAngle(dip); + setPhase(phase); -void Helix::setPhase(double val) { - mPhase = val; - mCosPhase = cos(mPhase); - mSinPhase = sin(mPhase); - if (fabs(mPhase) > M_PI) - mPhase = atan2(mSinPhase, mCosPhase); // force range [-pi,pi] + // + // Check for singularity and correct for negative curvature. + // May change mH and mPhase. Must therefore be set last. + // + setCurvature(c); + + // + // For the case B=0, h is ill defined. In the following we + // always assume h = +1. Since phase = psi - h * pi/2 + // we have to correct the phase in case h = -1. + // This assumes that the user uses the same h for phase + // as the one he passed to the constructor. + // + if (mSingularity && mH == -1) { + mH = +1; + setPhase(mPhase-M_PI); + } } -void Helix::setDipAngle(double val) { - mDipAngle = val; - mCosDipAngle = cos(mDipAngle); - mSinDipAngle = sin(mDipAngle); +void Helix::setCurvature(double val) +{ + if (val < 0) { + mCurvature = -val; + mH = -mH; + setPhase(mPhase+M_PI); + } + else + mCurvature = val; + + if (fabs(mCurvature) <= std::numeric_limits::epsilon()) + mSingularity = true; // straight line + else + mSingularity = false; // curved } -double Helix::xcenter() const { - if (mSingularity) - return 0; - else - return mOrigin.x - mCosPhase / mCurvature; +void Helix::setPhase(double val) +{ + mPhase = val; + mCosPhase = cos(mPhase); + mSinPhase = sin(mPhase); + if (fabs(mPhase) > M_PI) + mPhase = atan2(mSinPhase, mCosPhase); // force range [-pi,pi] } -double Helix::ycenter() const { - if (mSingularity) - return 0; - else - return mOrigin.y - mSinPhase / mCurvature; +void Helix::setDipAngle(double val) +{ + mDipAngle = val; + mCosDipAngle = cos(mDipAngle); + mSinDipAngle = sin(mDipAngle); } -double Helix::fudgePathLength(const edm4hep::Vector3f& p) const { - double s; - double dx = p.x - mOrigin.x; - double dy = p.y - mOrigin.y; +double Helix::xcenter() const +{ + if (mSingularity) + return 0; + else + return mOrigin.x-mCosPhase/mCurvature; +} - if (mSingularity) { - s = (dy * mCosPhase - dx * mSinPhase) / mCosDipAngle; - } else { - s = atan2(dy * mCosPhase - dx * mSinPhase, 1 / mCurvature + dx * mCosPhase + dy * mSinPhase) / - (mH * mCurvature * mCosDipAngle); - } - return s; +double Helix::ycenter() const +{ + if (mSingularity) + return 0; + else + return mOrigin.y-mSinPhase/mCurvature; } -double Helix::distance(const edm4hep::Vector3f& p, bool scanPeriods) const { - return edm4hep::utils::magnitude(this->at(pathLength(p, scanPeriods)) - p); +double Helix::fudgePathLength(const edm4hep::Vector3f& p) const +{ + double s; + double dx = p.x-mOrigin.x; + double dy = p.y-mOrigin.y; + + if (mSingularity) { + s = (dy*mCosPhase - dx*mSinPhase)/mCosDipAngle; + } + else { + s = atan2(dy*mCosPhase - dx*mSinPhase, + 1/mCurvature + dx*mCosPhase+dy*mSinPhase)/ + (mH*mCurvature*mCosDipAngle); + } + return s; } -double Helix::pathLength(const edm4hep::Vector3f& p, bool scanPeriods) const { - // - // Returns the path length at the distance of closest - // approach between the helix and point p. - // For the case of B=0 (straight line) the path length - // can be calculated analytically. For B>0 there is - // unfortunately no easy solution to the problem. - // Here we use the Newton method to find the root of the - // referring equation. The 'fudgePathLength' serves - // as a starting value. - // - double s; - double dx = p.x - mOrigin.x; - double dy = p.y - mOrigin.y; - double dz = p.z - mOrigin.z; - - if (mSingularity) { - s = mCosDipAngle * (mCosPhase * dy - mSinPhase * dx) + mSinDipAngle * dz; - } else { // - const double MaxPrecisionNeeded = edm4eic::unit::um; - const int MaxIterations = 100; +double Helix::distance(const edm4hep::Vector3f& p, bool scanPeriods) const +{ + return edm4hep::utils::magnitude(this->at(pathLength(p,scanPeriods))-p); +} +double Helix::pathLength(const edm4hep::Vector3f& p, bool scanPeriods) const +{ // - // The math is taken from Maple with C(expr,optimized) and - // some hand-editing. It is not very nice but efficient. + // Returns the path length at the distance of closest + // approach between the helix and point p. + // For the case of B=0 (straight line) the path length + // can be calculated analytically. For B>0 there is + // unfortunately no easy solution to the problem. + // Here we use the Newton method to find the root of the + // referring equation. The 'fudgePathLength' serves + // as a starting value. // - double t34 = mCurvature * mCosDipAngle * mCosDipAngle; - double t41 = mSinDipAngle * mSinDipAngle; - double t6, t7, t11, t12, t19; + double s; + double dx = p.x-mOrigin.x; + double dy = p.y-mOrigin.y; + double dz = p.z-mOrigin.z; - // - // Get a first guess by using the dca in 2D. Since - // in some extreme cases we might be off by n periods - // we add (subtract) periods in case we get any closer. - // - s = fudgePathLength(p); - - if (scanPeriods) { - double ds = period(); - int j, jmin = 0; - double d, dmin = edm4hep::utils::magnitude(at(s) - p); - for (j = 1; j < MaxIterations; j++) { - if ((d = edm4hep::utils::magnitude(at(s + j * ds) - p)) < dmin) { - dmin = d; - jmin = j; - } else - break; - } - for (j = -1; -j < MaxIterations; j--) { - if ((d = edm4hep::utils::magnitude(at(s + j * ds) - p)) < dmin) { - dmin = d; - jmin = j; - } else - break; - } - if (jmin) - s += jmin * ds; + if (mSingularity) { + s = mCosDipAngle*(mCosPhase*dy-mSinPhase*dx) + + mSinDipAngle*dz; } - - // - // Newtons method: - // Stops after MaxIterations iterations or if the required - // precision is obtained. Whatever comes first. - // - double sOld = s; - for (int i = 0; i < MaxIterations; i++) { - t6 = mPhase + s * mH * mCurvature * mCosDipAngle; - t7 = cos(t6); - t11 = dx - (1 / mCurvature) * (t7 - mCosPhase); - t12 = sin(t6); - t19 = dy - (1 / mCurvature) * (t12 - mSinPhase); - s -= (t11 * t12 * mH * mCosDipAngle - t19 * t7 * mH * mCosDipAngle - - (dz - s * mSinDipAngle) * mSinDipAngle) / - (t12 * t12 * mCosDipAngle * mCosDipAngle + t11 * t7 * t34 + - t7 * t7 * mCosDipAngle * mCosDipAngle + t19 * t12 * t34 + t41); - if (fabs(sOld - s) < MaxPrecisionNeeded) - break; - sOld = s; + else { // + const double MaxPrecisionNeeded = edm4eic::unit::um; + const int MaxIterations = 100; + + // + // The math is taken from Maple with C(expr,optimized) and + // some hand-editing. It is not very nice but efficient. + // + double t34 = mCurvature*mCosDipAngle*mCosDipAngle; + double t41 = mSinDipAngle*mSinDipAngle; + double t6, t7, t11, t12, t19; + + // + // Get a first guess by using the dca in 2D. Since + // in some extreme cases we might be off by n periods + // we add (subtract) periods in case we get any closer. + // + s = fudgePathLength(p); + + if (scanPeriods) { + double ds = period(); + int j, jmin = 0; + double d, dmin = edm4hep::utils::magnitude(at(s) - p); + for(j=1; j::max(); - else - return fabs(2 * M_PI / (mH * mCurvature * mCosDipAngle)); +double Helix::period() const +{ + if (mSingularity) + return std::numeric_limits::max(); + else + return fabs(2*M_PI/(mH*mCurvature*mCosDipAngle)); } -std::pair Helix::pathLength(double r) const { - std::pair value; - std::pair VALUE(999999999., 999999999.); - // - // The math is taken from Maple with C(expr,optimized) and - // some hand-editing. It is not very nice but efficient. - // 'first' is the smallest of the two solutions (may be negative) - // 'second' is the other. - // - if (mSingularity) { - double t1 = mCosDipAngle * (mOrigin.x * mSinPhase - mOrigin.y * mCosPhase); - double t12 = mOrigin.y * mOrigin.y; - double t13 = mCosPhase * mCosPhase; - double t15 = r * r; - double t16 = mOrigin.x * mOrigin.x; - double t20 = - -mCosDipAngle * mCosDipAngle * - (2.0 * mOrigin.x * mSinPhase * mOrigin.y * mCosPhase + t12 - t12 * t13 - t15 + t13 * t16); - if (t20 < 0.) - return VALUE; - t20 = ::sqrt(t20); - value.first = (t1 - t20) / (mCosDipAngle * mCosDipAngle); - value.second = (t1 + t20) / (mCosDipAngle * mCosDipAngle); - } else { - double t1 = mOrigin.y * mCurvature; - double t2 = mSinPhase; - double t3 = mCurvature * mCurvature; - double t4 = mOrigin.y * t2; - double t5 = mCosPhase; - double t6 = mOrigin.x * t5; - double t8 = mOrigin.x * mOrigin.x; - double t11 = mOrigin.y * mOrigin.y; - double t14 = r * r; - double t15 = t14 * mCurvature; - double t17 = t8 * t8; - double t19 = t11 * t11; - double t21 = t11 * t3; - double t23 = t5 * t5; - double t32 = t14 * t14; - double t35 = t14 * t3; - double t38 = 8.0 * t4 * t6 - 4.0 * t1 * t2 * t8 - 4.0 * t11 * mCurvature * t6 + 4.0 * t15 * t6 + - t17 * t3 + t19 * t3 + 2.0 * t21 * t8 + 4.0 * t8 * t23 - - 4.0 * t8 * mOrigin.x * mCurvature * t5 - 4.0 * t11 * t23 - - 4.0 * t11 * mOrigin.y * mCurvature * t2 + 4.0 * t11 - 4.0 * t14 + t32 * t3 + - 4.0 * t15 * t4 - 2.0 * t35 * t11 - 2.0 * t35 * t8; - double t40 = (-t3 * t38); - if (t40 < 0.) - return VALUE; - t40 = ::sqrt(t40); - - double t43 = mOrigin.x * mCurvature; - double t45 = 2.0 * t5 - t35 + t21 + 2.0 - 2.0 * t1 * t2 - 2.0 * t43 - 2.0 * t43 * t5 + t8 * t3; - double t46 = mH * mCosDipAngle * mCurvature; - - value.first = (-mPhase + 2.0 * atan((-2.0 * t1 + 2.0 * t2 + t40) / t45)) / t46; - value.second = -(mPhase + 2.0 * atan((2.0 * t1 - 2.0 * t2 + t40) / t45)) / t46; - +std::pair Helix::pathLength(double r) const +{ + std::pair value; + std::pair VALUE(999999999.,999999999.); // - // Solution can be off by +/- one period, select smallest + // The math is taken from Maple with C(expr,optimized) and + // some hand-editing. It is not very nice but efficient. + // 'first' is the smallest of the two solutions (may be negative) + // 'second' is the other. // - double p = period(); - if (!std::isnan(value.first)) { - if (fabs(value.first - p) < fabs(value.first)) - value.first = value.first - p; - else if (fabs(value.first + p) < fabs(value.first)) - value.first = value.first + p; + if (mSingularity) { + double t1 = mCosDipAngle*(mOrigin.x*mSinPhase-mOrigin.y*mCosPhase); + double t12 = mOrigin.y*mOrigin.y; + double t13 = mCosPhase*mCosPhase; + double t15 = r*r; + double t16 = mOrigin.x*mOrigin.x; + double t20 = -mCosDipAngle*mCosDipAngle*(2.0*mOrigin.x*mSinPhase*mOrigin.y*mCosPhase + + t12-t12*t13-t15+t13*t16); + if (t20<0.) return VALUE; + t20 = ::sqrt(t20); + value.first = (t1-t20)/(mCosDipAngle*mCosDipAngle); + value.second = (t1+t20)/(mCosDipAngle*mCosDipAngle); } - if (!std::isnan(value.second)) { - if (fabs(value.second - p) < fabs(value.second)) - value.second = value.second - p; - else if (fabs(value.second + p) < fabs(value.second)) - value.second = value.second + p; + else { + double t1 = mOrigin.y*mCurvature; + double t2 = mSinPhase; + double t3 = mCurvature*mCurvature; + double t4 = mOrigin.y*t2; + double t5 = mCosPhase; + double t6 = mOrigin.x*t5; + double t8 = mOrigin.x*mOrigin.x; + double t11 = mOrigin.y*mOrigin.y; + double t14 = r*r; + double t15 = t14*mCurvature; + double t17 = t8*t8; + double t19 = t11*t11; + double t21 = t11*t3; + double t23 = t5*t5; + double t32 = t14*t14; + double t35 = t14*t3; + double t38 = 8.0*t4*t6 - 4.0*t1*t2*t8 - 4.0*t11*mCurvature*t6 + + 4.0*t15*t6 + t17*t3 + t19*t3 + 2.0*t21*t8 + 4.0*t8*t23 - + 4.0*t8*mOrigin.x*mCurvature*t5 - 4.0*t11*t23 - + 4.0*t11*mOrigin.y*mCurvature*t2 + 4.0*t11 - 4.0*t14 + + t32*t3 + 4.0*t15*t4 - 2.0*t35*t11 - 2.0*t35*t8; + double t40 = (-t3*t38); + if (t40<0.) return VALUE; + t40 = ::sqrt(t40); + + double t43 = mOrigin.x*mCurvature; + double t45 = 2.0*t5 - t35 + t21 + 2.0 - 2.0*t1*t2 -2.0*t43 - 2.0*t43*t5 + t8*t3; + double t46 = mH*mCosDipAngle*mCurvature; + + value.first = (-mPhase + 2.0*atan((-2.0*t1 + 2.0*t2 + t40)/t45))/t46; + value.second = -(mPhase + 2.0*atan((2.0*t1 - 2.0*t2 + t40)/t45))/t46; + + // + // Solution can be off by +/- one period, select smallest + // + double p = period(); + if (! std::isnan(value.first)) { + if (fabs(value.first-p) < fabs(value.first)) value.first = value.first-p; + else if (fabs(value.first+p) < fabs(value.first)) value.first = value.first+p; + } + if (! std::isnan(value.second)) { + if (fabs(value.second-p) < fabs(value.second)) value.second = value.second-p; + else if (fabs(value.second+p) < fabs(value.second)) value.second = value.second+p; + } } - } - if (value.first > value.second) - std::swap(value.first, value.second); - return (value); -} - -std::pair Helix::pathLength(double r, double x, double y) { - double x0 = mOrigin.x; - double y0 = mOrigin.y; - mOrigin.x = x0 - x; - mOrigin.y = y0 - y; - std::pair result = this->pathLength(r); - mOrigin.x = x0; - mOrigin.y = y0; - return result; + if (value.first > value.second) + std::swap(value.first,value.second); + return(value); } -double Helix::pathLength(const edm4hep::Vector3f& r, const edm4hep::Vector3f& n) const { - // - // Vector 'r' defines the position of the center and - // vector 'n' the normal vector of the plane. - // For a straight line there is a simple analytical - // solution. For curvatures > 0 the root is determined - // by Newton method. In case no valid s can be found - // the max. largest value for s is returned. - // - double s; - - if (mSingularity) { - double t = n.z * mSinDipAngle + n.y * mCosDipAngle * mCosPhase - n.x * mCosDipAngle * mSinPhase; - if (t == 0) - s = NoSolution; - else - s = ((r - mOrigin) * n) / t; - } else { - const double MaxPrecisionNeeded = edm4eic::unit::um; - const int MaxIterations = 20; - - double A = mCurvature * ((mOrigin - r) * n) - n.x * mCosPhase - n.y * mSinPhase; - double t = mH * mCurvature * mCosDipAngle; - double u = n.z * mCurvature * mSinDipAngle; - - double a, f, fp; - double sOld = s = 0; - double shiftOld = 0; - double shift; - // (cos(angMax)-1)/angMax = 0.1 - const double angMax = 0.21; - double deltas = fabs(angMax / (mCurvature * mCosDipAngle)); - // dampingFactor = exp(-0.5); - // double dampingFactor = 0.60653; - int i; - - for (i = 0; i < MaxIterations; i++) { - a = t * s + mPhase; - double sina = sin(a); - double cosa = cos(a); - f = A + n.x * cosa + n.y * sina + u * s; - fp = -n.x * sina * t + n.y * cosa * t + u; - if (fabs(fp) * deltas <= fabs(f)) { //too big step - int sgn = 1; - if (fp < 0.) - sgn = -sgn; - if (f < 0.) - sgn = -sgn; - shift = sgn * deltas; - if (shift < 0) - shift *= 0.9; // don't get stuck shifting +/-deltas - } else { - shift = f / fp; - } - s -= shift; - shiftOld = shift; - if (fabs(sOld - s) < MaxPrecisionNeeded) - break; - sOld = s; - } - if (i == MaxIterations) - return NoSolution; - } - return s; +std::pair Helix::pathLength(double r, double x, double y) +{ + double x0 = mOrigin.x; + double y0 = mOrigin.y; + mOrigin.x = x0-x; + mOrigin.y = y0-y; + std::pair result = this->pathLength(r); + mOrigin.x = x0; + mOrigin.y = y0; + return result; } -std::pair Helix::pathLengths(const Helix& h, double minStepSize, - double minRange) const { - // - // Cannot handle case where one is a helix - // and the other one is a straight line. - // - if (mSingularity != h.mSingularity) - return std::pair(NoSolution, NoSolution); - - double s1, s2; - - if (mSingularity) { - // - // Analytic solution - // - edm4hep::Vector3f dv = h.mOrigin - mOrigin; - edm4hep::Vector3f a(-mCosDipAngle * mSinPhase, mCosDipAngle * mCosPhase, mSinDipAngle); - edm4hep::Vector3f b(-h.mCosDipAngle * h.mSinPhase, h.mCosDipAngle * h.mCosPhase, - h.mSinDipAngle); - double ab = a * b; - double g = dv * a; - double k = dv * b; - s2 = (k - ab * g) / (ab * ab - 1.); - s1 = g + s2 * ab; - return std::pair(s1, s2); - } else { +double Helix::pathLength(const edm4hep::Vector3f& r, + const edm4hep::Vector3f& n) const +{ // - // First step: get dca in the xy-plane as start value + // Vector 'r' defines the position of the center and + // vector 'n' the normal vector of the plane. + // For a straight line there is a simple analytical + // solution. For curvatures > 0 the root is determined + // by Newton method. In case no valid s can be found + // the max. largest value for s is returned. // - double dx = h.xcenter() - xcenter(); - double dy = h.ycenter() - ycenter(); - double dd = ::sqrt(dx * dx + dy * dy); - double r1 = 1 / curvature(); - double r2 = 1 / h.curvature(); - - double cosAlpha = (r1 * r1 + dd * dd - r2 * r2) / (2 * r1 * dd); - double s; - double x, y; - if (fabs(cosAlpha) < 1) { // two solutions - double sinAlpha = sin(acos(cosAlpha)); - x = xcenter() + r1 * (cosAlpha * dx - sinAlpha * dy) / dd; - y = ycenter() + r1 * (sinAlpha * dx + cosAlpha * dy) / dd; - s = pathLength(x, y); - x = xcenter() + r1 * (cosAlpha * dx + sinAlpha * dy) / dd; - y = ycenter() + r1 * (cosAlpha * dy - sinAlpha * dx) / dd; - double a = pathLength(x, y); - if (h.distance(at(a)) < h.distance(at(s))) - s = a; - } else { // no intersection (or exactly one) - int rsign = ((r2 - r1) > dd ? -1 : 1); // set -1 when *this* helix is - // completely contained in the other - x = xcenter() + rsign * r1 * dx / dd; - y = ycenter() + rsign * r1 * dy / dd; - s = pathLength(x, y); - } - // - // Second step: scan in decreasing intervals around seed 's' - // minRange and minStepSize are passed as arguments to the method. - // They have default values defined in the header file. - // - double dmin = h.distance(at(s)); - double range = std::max(2 * dmin, minRange); - double ds = range / 10; - double slast = -999999, ss, d; - s1 = s - range / 2.; - s2 = s + range / 2.; - - while (ds > minStepSize) { - for (ss = s1; ss < s2 + ds; ss += ds) { - d = h.distance(at(ss)); - if (d < dmin) { - dmin = d; - s = ss; - } - slast = ss; - } - // - // In the rare cases where the minimum is at the - // the border of the current range we shift the range - // and start all over, i.e we do not decrease 'ds'. - // Else we decrease the search intervall around the - // current minimum and redo the scan in smaller steps. - // - if (s == s1) { - d = 0.8 * (s2 - s1); - s1 -= d; - s2 -= d; - } else if (s == slast) { - d = 0.8 * (s2 - s1); - s1 += d; - s2 += d; - } else { - s1 = s - ds; - s2 = s + ds; - ds /= 10; - } + if (mSingularity) { + double t = n.z*mSinDipAngle + + n.y*mCosDipAngle*mCosPhase - + n.x*mCosDipAngle*mSinPhase; + if (t == 0) + s = NoSolution; + else + s = ((r - mOrigin)*n)/t; } - return std::pair(s, h.pathLength(at(s))); - } + else { + const double MaxPrecisionNeeded = edm4eic::unit::um; + const int MaxIterations = 20; + + double A = mCurvature*((mOrigin - r)*n) - + n.x*mCosPhase - + n.y*mSinPhase; + double t = mH*mCurvature*mCosDipAngle; + double u = n.z*mCurvature*mSinDipAngle; + + double a, f, fp; + double sOld = s = 0; + double shiftOld = 0; + double shift; +// (cos(angMax)-1)/angMax = 0.1 + const double angMax = 0.21; + double deltas = fabs(angMax/(mCurvature*mCosDipAngle)); +// dampingFactor = exp(-0.5); +// double dampingFactor = 0.60653; + int i; + + for (i=0; i Helix::pathLengths(const Helix& h, double minStepSize, - double minRange) const { - // - // Cannot handle case where one is a helix - // and the other one is a straight line. - // - if (mSingularity != h.mSingularity) - return std::pair(NoSolution, NoSolution); - - double s1, s2; - - if (mSingularity) { - // - // Analytic solution - // - edm4hep::Vector3f dv = h.mOrigin - mOrigin; - edm4hep::Vector3f a(-mCosDipAngle * mSinPhase, mCosDipAngle * mCosPhase, mSinDipAngle); - edm4hep::Vector3f b(-h.mCosDipAngle * h.mSinPhase, h.mCosDipAngle * h.mCosPhase, - h.mSinDipAngle); - double ab = a * b; - double g = dv * a; - double k = dv * b; - s2 = (k - ab * g) / (ab * ab - 1.); - s1 = g + s2 * ab; - return std::pair(s1, s2); - } else { +std::pair +Helix::pathLengths(const Helix& h, double minStepSize, double minRange) const +{ // - // First step: get dca in the xy-plane as start value + // Cannot handle case where one is a helix + // and the other one is a straight line. // - double dx = h.xcenter() - xcenter(); - double dy = h.ycenter() - ycenter(); - double dd = ::sqrt(dx * dx + dy * dy); - double r1 = 1 / curvature(); - double r2 = 1 / h.curvature(); - - double cosAlpha = (r1 * r1 + dd * dd - r2 * r2) / (2 * r1 * dd); - - double s; - double x, y; - if (fabs(cosAlpha) < 1) { // two solutions - double sinAlpha = sin(acos(cosAlpha)); - x = xcenter() + r1 * (cosAlpha * dx - sinAlpha * dy) / dd; - y = ycenter() + r1 * (sinAlpha * dx + cosAlpha * dy) / dd; - s = pathLength(x, y); - x = xcenter() + r1 * (cosAlpha * dx + sinAlpha * dy) / dd; - y = ycenter() + r1 * (cosAlpha * dy - sinAlpha * dx) / dd; - double a = pathLength(x, y); - if (h.distance(at(a)) < h.distance(at(s))) - s = a; - } else { // no intersection (or exactly one) - int rsign = ((r2 - r1) > dd ? -1 : 1); // set -1 when *this* helix is - // completely contained in the other - x = xcenter() + rsign * r1 * dx / dd; - y = ycenter() + rsign * r1 * dy / dd; - s = pathLength(x, y); + if (mSingularity != h.mSingularity) + return std::pair(NoSolution, NoSolution); + + double s1, s2; + + if (mSingularity) { + // + // Analytic solution + // + edm4hep::Vector3f dv = h.mOrigin - mOrigin; + edm4hep::Vector3f a(-mCosDipAngle*mSinPhase, + mCosDipAngle*mCosPhase, + mSinDipAngle); + edm4hep::Vector3f b(-h.mCosDipAngle*h.mSinPhase, + h.mCosDipAngle*h.mCosPhase, + h.mSinDipAngle); + double ab = a*b; + double g = dv*a; + double k = dv*b; + s2 = (k-ab*g)/(ab*ab-1.); + s1 = g+s2*ab; + return std::pair(s1, s2); } - - // - // Second step: scan in decreasing intervals around seed 's' - // minRange and minStepSize are passed as arguments to the method. - // They have default values defined in the header file. - // - double dmin = h.distance(at(s)); - double range = std::max(2 * dmin, minRange); - double ds = range / 10; - double slast = -999999, ss, d; - s1 = s - range / 2.; - s2 = s + range / 2.; - - while (ds > minStepSize) { - for (ss = s1; ss < s2 + ds; ss += ds) { - d = h.distance(at(ss)); - if (d < dmin) { - dmin = d; - s = ss; + else { + // + // First step: get dca in the xy-plane as start value + // + double dx = h.xcenter() - xcenter(); + double dy = h.ycenter() - ycenter(); + double dd = ::sqrt(dx*dx + dy*dy); + double r1 = 1/curvature(); + double r2 = 1/h.curvature(); + + double cosAlpha = (r1*r1 + dd*dd - r2*r2)/(2*r1*dd); + + double s; + double x, y; + if (fabs(cosAlpha) < 1) { // two solutions + double sinAlpha = sin(acos(cosAlpha)); + x = xcenter() + r1*(cosAlpha*dx - sinAlpha*dy)/dd; + y = ycenter() + r1*(sinAlpha*dx + cosAlpha*dy)/dd; + s = pathLength(x, y); + x = xcenter() + r1*(cosAlpha*dx + sinAlpha*dy)/dd; + y = ycenter() + r1*(cosAlpha*dy - sinAlpha*dx)/dd; + double a = pathLength(x, y); + if (h.distance(at(a)) < h.distance(at(s))) s = a; + } + else { // no intersection (or exactly one) + int rsign = ((r2-r1) > dd ? -1 : 1); // set -1 when *this* helix is + // completely contained in the other + x = xcenter() + rsign*r1*dx/dd; + y = ycenter() + rsign*r1*dy/dd; + s = pathLength(x, y); } - slast = ss; - } - // - // In the rare cases where the minimum is at the - // the border of the current range we shift the range - // and start all over, i.e we do not decrease 'ds'. - // Else we decrease the search interval around the - // current minimum and redo the scan in smaller steps. - // - if (s == s1) { - d = 0.8 * (s2 - s1); - s1 -= d; - s2 -= d; - } else if (s == slast) { - d = 0.8 * (s2 - s1); - s1 += d; - s2 += d; - } else { - s1 = s - ds; - s2 = s + ds; - ds /= 10; - } + + // + // Second step: scan in decreasing intervals around seed 's' + // minRange and minStepSize are passed as arguments to the method. + // They have default values defined in the header file. + // + double dmin = h.distance(at(s)); + double range = std::max(2*dmin, minRange); + double ds = range/10; + double slast=-999999, ss, d; + s1 = s - range/2.; + s2 = s + range/2.; + + while (ds > minStepSize) { + for (ss=s1; ss(s, h.pathLength(at(s))); } - return std::pair(s, h.pathLength(at(s))); - } } -void Helix::moveOrigin(double s) { - if (mSingularity) - mOrigin = at(s); - else { - edm4hep::Vector3f newOrigin = at(s); - double newPhase = atan2(newOrigin.y - ycenter(), newOrigin.x - xcenter()); - mOrigin = newOrigin; - setPhase(newPhase); - } ->>>>>>> 5c2ad0d0 (fixed typos in Helix function (in commenting area)) + +void Helix::moveOrigin(double s) +{ + if (mSingularity) + mOrigin = at(s); + else { + edm4hep::Vector3f newOrigin = at(s); + double newPhase = atan2(newOrigin.y - ycenter(), + newOrigin.x - xcenter()); + mOrigin = newOrigin; + setPhase(newPhase); + } } /* int operator== (const Helix& a, const Helix& b) @@ -646,88 +560,100 @@ int operator== (const Helix& a, const Helix& b) // Checks for numerical identity only ! // return (a.origin() == b.origin() && - a.dipAngle() == b.dipAngle() && - a.curvature() == b.curvature() && - a.phase() == b.phase() && - a.h() == b.h()); + a.dipAngle() == b.dipAngle() && + a.curvature() == b.curvature() && + a.phase() == b.phase() && + a.h() == b.h()); } int operator!= (const Helix& a, const Helix& b) {return !(a == b);} */ //std::ostream& operator<<(std::ostream& os, const Helix& h) -void Helix::Print() const { - std::cout << "(" - << "curvature = " << mCurvature << ", " - << "dip angle = " << mDipAngle << ", " - << "phase = " << mPhase << ", " - << "h = " << mH << ", " - << "origin = " << mOrigin.x << " " << mOrigin.y << " " << mOrigin.z << ")" << std::endl; +void Helix::Print() const +{ + std::cout << "(" + << "curvature = " << mCurvature << ", " + << "dip angle = " << mDipAngle << ", " + << "phase = " << mPhase << ", " + << "h = " << mH << ", " + << "origin = " << mOrigin.x << " " << mOrigin.y << " " << mOrigin.z << ")" + << std::endl; } -edm4hep::Vector3f Helix::momentum(double B) const { - if (mSingularity) - return (edm4hep::Vector3f(0, 0, 0)); - else { - double pt = edm4eic::unit::GeV * - fabs(dd4hep::c_light * dd4hep::nanosecond / dd4hep::meter * B / dd4hep::tesla) / - (fabs(mCurvature) * dd4hep::meter); - - return (edm4hep::Vector3f(pt * cos(mPhase + mH * M_PI / 2), // pos part pos field - pt * sin(mPhase + mH * M_PI / 2), pt * tan(mDipAngle))); - } +edm4hep::Vector3f Helix::momentum(double B) const +{ + if (mSingularity) + return(edm4hep::Vector3f(0,0,0)); + else { + double pt = edm4eic::unit::GeV*fabs(dd4hep::c_light*dd4hep::nanosecond/dd4hep::meter*B/dd4hep::tesla)/(fabs(mCurvature)*dd4hep::meter); + + return (edm4hep::Vector3f(pt*cos(mPhase+mH*M_PI/2), // pos part pos field + pt*sin(mPhase+mH*M_PI/2), + pt*tan(mDipAngle))); + } } -edm4hep::Vector3f Helix::momentumAt(double S, double B) const { - // Obtain phase-shifted momentum from phase-shift of origin - Helix tmp(*this); - tmp.moveOrigin(S); - return tmp.momentum(B); +edm4hep::Vector3f Helix::momentumAt(double S, double B) const +{ + // Obtain phase-shifted momentum from phase-shift of origin + Helix tmp(*this); + tmp.moveOrigin(S); + return tmp.momentum(B); } -int Helix::charge(double B) const { return (B > 0 ? -mH : mH); } - -double Helix::geometricSignedDistance(double x, double y) { - // Geometric signed distance - double thePath = this->pathLength(x, y); - edm4hep::Vector3f DCA2dPosition = this->at(thePath); - DCA2dPosition.z = 0; - edm4hep::Vector3f position(x, y, 0); - edm4hep::Vector3f DCAVec = (DCA2dPosition - position); - edm4hep::Vector3f momVec; - // Deal with straight tracks - if (this->mSingularity) { - momVec = this->at(1) - this->at(0); - momVec.z = 0; - } else { - momVec = this->momentumAt( - thePath, 1. / dd4hep::tesla); // Don't care about Bmag. Helicity is what matters. - momVec.z = 0; - } +int Helix::charge(double B) const +{ + return (B > 0 ? -mH : mH); +} - double cross = DCAVec.x * momVec.y - DCAVec.y * momVec.x; - double theSign = (cross >= 0) ? 1. : -1.; - return theSign * edm4hep::utils::magnitudeTransverse(DCAVec); +double Helix::geometricSignedDistance(double x, double y) +{ + // Geometric signed distance + double thePath = this->pathLength(x,y); + edm4hep::Vector3f DCA2dPosition = this->at(thePath); + DCA2dPosition.z = 0; + edm4hep::Vector3f position(x,y,0); + edm4hep::Vector3f DCAVec = (DCA2dPosition-position); + edm4hep::Vector3f momVec; + // Deal with straight tracks + if (this->mSingularity) { + momVec = this->at(1)- this->at(0); + momVec.z = 0; + } + else { + momVec = this->momentumAt(thePath,1./dd4hep::tesla); // Don't care about Bmag. Helicity is what matters. + momVec.z = 0; + } + + double cross = DCAVec.x*momVec.y - DCAVec.y*momVec.x; + double theSign = (cross>=0) ? 1. : -1.; + return theSign*edm4hep::utils::magnitudeTransverse(DCAVec); } -double Helix::curvatureSignedDistance(double x, double y) { - // Protect against mH = 0 or zero field - if (this->mSingularity || abs(this->mH) <= 0) { - return (this->geometricSignedDistance(x, y)); - } else { - return (this->geometricSignedDistance(x, y)) / (this->mH); - } +double Helix::curvatureSignedDistance(double x, double y) +{ + // Protect against mH = 0 or zero field + if (this->mSingularity || abs(this->mH)<=0) { + return (this->geometricSignedDistance(x,y)); + } + else { + return (this->geometricSignedDistance(x,y))/(this->mH); + } + } -double Helix::geometricSignedDistance(const edm4hep::Vector3f& pos) { - double sdca2d = this->geometricSignedDistance(pos.x, pos.y); - double theSign = (sdca2d >= 0) ? 1. : -1.; - return (this->distance(pos)) * theSign; +double Helix::geometricSignedDistance(const edm4hep::Vector3f& pos) +{ + double sdca2d = this->geometricSignedDistance(pos.x,pos.y); + double theSign = (sdca2d>=0) ? 1. : -1.; + return (this->distance(pos))*theSign; } -double Helix::curvatureSignedDistance(const edm4hep::Vector3f& pos) { - double sdca2d = this->curvatureSignedDistance(pos.x, pos.y); - double theSign = (sdca2d >= 0) ? 1. : -1.; - return (this->distance(pos)) * theSign; +double Helix::curvatureSignedDistance(const edm4hep::Vector3f& pos) +{ + double sdca2d = this->curvatureSignedDistance(pos.x,pos.y); + double theSign = (sdca2d>=0) ? 1. : -1.; + return (this->distance(pos))*theSign; } -} // namespace eicrecon +} // namespace eicrecon \ No newline at end of file diff --git a/src/algorithms/reco/Helix.h b/src/algorithms/reco/Helix.h index 2f30c1d6ce..437336343d 100644 --- a/src/algorithms/reco/Helix.h +++ b/src/algorithms/reco/Helix.h @@ -21,136 +21,119 @@ namespace eicrecon { class Helix { - bool mSingularity; // true for straight line case (B=0) - edm4hep::Vector3f mOrigin; - double mDipAngle; - double mCurvature; - double mPhase; - int mH; // -sign(q*B); - - double mCosDipAngle; - double mSinDipAngle; - double mCosPhase; - double mSinPhase; + bool mSingularity; // true for straight line case (B=0) + edm4hep::Vector3f mOrigin; + double mDipAngle; + double mCurvature; + double mPhase; + int mH; // -sign(q*B); + + double mCosDipAngle; + double mSinDipAngle; + double mCosPhase; + double mSinPhase; public: - /// curvature, dip angle, phase, origin, h - Helix(const double c, const double dip, const double phase, const edm4hep::Vector3f& o, - const int h = -1); - - /// momentum, origin, b_field, charge - Helix(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q); - - /// ReconstructParticle, b field - Helix(const edm4eic::ReconstructedParticle& p, const double b_field); - - ~Helix() = default; - - double dipAngle() const; - double curvature() const; /// 1/R in xy-plane - double phase() const; /// aziumth in xy-plane measured from ring center - double xcenter() const; /// x-center of circle in xy-plane - double ycenter() const; /// y-center of circle in xy-plane - int h() const; /// -sign(q*B); - - const edm4hep::Vector3f& origin() const; /// starting point - - void setParameters(double c, double dip, double phase, const edm4hep::Vector3f& o, int h); - - void setParameters(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, - const int q); - - /// edm4eic::TrackParameters, b field - void setParameters(const edm4eic::TrackParameters& trk, const double b_field); - - /// coordinates of helix at point s - double x(double s) const; - double y(double s) const; - double z(double s) const; - - edm4hep::Vector3f at(double s) const; - - /// pointing vector of helix at point s - double cx(double s) const; - double cy(double s) const; - double cz(double s = 0) const; - - edm4hep::Vector3f cat(double s) const; - -<<<<<<< HEAD - /// returns period length of helix - double period() const; - - /// path length at given r (cylindrical r) - std::pair pathLength(double r) const; - - /// path length at given r (cylindrical r, cylinder axis at x,y) - std::pair pathLength(double r, double x, double y); - - /// path length at distance of closest approach to a given point - double pathLength(const edm4hep::Vector3f& p, bool scanPeriods = true) const; - - /// path length at intersection with plane - double pathLength(const edm4hep::Vector3f& r, const edm4hep::Vector3f& n) const; - - /// path length at distance of closest approach in the xy-plane to a given point - double pathLength(double x, double y) const; - - /// path lengths at dca between two helices - std::pair pathLengths(const Helix&, double minStepSize = 10 * edm4eic::unit::um, - double minRange = 10 * edm4eic::unit::cm) const; - - /// minimal distance between point and helix - double distance(const edm4hep::Vector3f& p, bool scanPeriods = true) const; - - /// checks for valid parametrization - bool valid(double world = 1.e+5) const { return !bad(world); } - int bad(double world = 1.e+5) const; - - /// move the origin along the helix to s which becomes then s=0 - virtual void moveOrigin(double s); - - static const double NoSolution; - - void setCurvature(double); /// performs also various checks - void setPhase(double); - void setDipAngle(double); - - /// value of S where distance in x-y plane is minimal - double fudgePathLength(const edm4hep::Vector3f&) const; - - // Requires: signed Magnetic Field - edm4hep::Vector3f momentum(double) const; // returns the momentum at origin - edm4hep::Vector3f momentumAt(double, double) const; // returns momemtum at S - int charge(double) const; // returns charge of particle - // 2d DCA to x,y point signed relative to curvature - double curvatureSignedDistance(double x, double y); - // 2d DCA to x,y point signed relative to rotation - double geometricSignedDistance(double x, double y); - // 3d DCA to 3d point signed relative to curvature - double curvatureSignedDistance(const edm4hep::Vector3f&); - // 3d DCA to 3d point signed relative to rotation - double geometricSignedDistance(const edm4hep::Vector3f&); - - // - void Print() const; -======= - // Requires: signed Magnetic Field - edm4hep::Vector3f momentum(double) const; // returns the momentum at origin - edm4hep::Vector3f momentumAt(double, double) const; // returns momentum at S - int charge(double) const; // returns charge of particle - // 2d DCA to x,y point signed relative to curvature - double curvatureSignedDistance(double x, double y); - // 2d DCA to x,y point signed relative to rotation - double geometricSignedDistance(double x, double y); - // 3d DCA to 3d point signed relative to curvature - double curvatureSignedDistance(const edm4hep::Vector3f&); - // 3d DCA to 3d point signed relative to rotation - double geometricSignedDistance(const edm4hep::Vector3f&); - - // - void Print() const; ->>>>>>> 5c2ad0d0 (fixed typos in Helix function (in commenting area)) + /// curvature, dip angle, phase, origin, h + Helix(const double c, const double dip, const double phase, const edm4hep::Vector3f& o, const int h=-1); + + /// momentum, origin, b_field, charge + Helix(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q); + + /// ReconstructParticle, b field + Helix(const edm4eic::ReconstructedParticle& p, const double b_field); + + ~Helix() = default; + + double dipAngle() const; + double curvature() const; /// 1/R in xy-plane + double phase() const; /// aziumth in xy-plane measured from ring center + double xcenter() const; /// x-center of circle in xy-plane + double ycenter() const; /// y-center of circle in xy-plane + int h() const; /// -sign(q*B); + + const edm4hep::Vector3f& origin() const; /// starting point + + void setParameters(double c, double dip, double phase, const edm4hep::Vector3f& o, int h); + + void setParameters(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q); + + /// edm4eic::TrackParameters, b field + void setParameters(const edm4eic::TrackParameters& trk, const double b_field); + + /// coordinates of helix at point s + double x(double s) const; + double y(double s) const; + double z(double s) const; + + edm4hep::Vector3f at(double s) const; + + /// pointing vector of helix at point s + double cx(double s) const; + double cy(double s) const; + double cz(double s = 0) const; + + edm4hep::Vector3f cat(double s) const; + + /// returns period length of helix + double period() const; + + /// path length at given r (cylindrical r) + std::pair pathLength(double r) const; + + /// path length at given r (cylindrical r, cylinder axis at x,y) + std::pair pathLength(double r, double x, double y); + + /// path length at distance of closest approach to a given point + double pathLength(const edm4hep::Vector3f& p, bool scanPeriods = true) const; + + /// path length at intersection with plane + double pathLength(const edm4hep::Vector3f& r, + const edm4hep::Vector3f& n) const; + + /// path length at distance of closest approach in the xy-plane to a given point + double pathLength(double x, double y) const; + + /// path lengths at dca between two helices + std::pair pathLengths(const Helix&, + double minStepSize = 10*edm4eic::unit::um, + double minRange = 10*edm4eic::unit::cm) const; + + /// minimal distance between point and helix + double distance(const edm4hep::Vector3f& p, bool scanPeriods = true) const; + + /// checks for valid parametrization + bool valid(double world = 1.e+5) const {return !bad(world);} + int bad(double world = 1.e+5) const; + + /// move the origin along the helix to s which becomes then s=0 + virtual void moveOrigin(double s); + + static const double NoSolution; + + + void setCurvature(double); /// performs also various checks + void setPhase(double); + void setDipAngle(double); + + /// value of S where distance in x-y plane is minimal + double fudgePathLength(const edm4hep::Vector3f&) const; + + // Requires: signed Magnetic Field + edm4hep::Vector3f momentum(double) const; // returns the momentum at origin + edm4hep::Vector3f momentumAt(double, double) const; // returns momentum at S + int charge(double) const; // returns charge of particle + // 2d DCA to x,y point signed relative to curvature + double curvatureSignedDistance(double x, double y) ; + // 2d DCA to x,y point signed relative to rotation + double geometricSignedDistance(double x, double y) ; + // 3d DCA to 3d point signed relative to curvature + double curvatureSignedDistance(const edm4hep::Vector3f&) ; + // 3d DCA to 3d point signed relative to rotation + double geometricSignedDistance(const edm4hep::Vector3f&) ; + + // + void Print() const; }; // end class Helix @@ -164,83 +147,92 @@ class Helix { // // Inline functions // -inline int Helix::h() const { return mH; } +inline int Helix::h() const {return mH;} -inline double Helix::dipAngle() const { return mDipAngle; } +inline double Helix::dipAngle() const {return mDipAngle;} -inline double Helix::curvature() const { return mCurvature; } +inline double Helix::curvature() const {return mCurvature;} -inline double Helix::phase() const { return mPhase; } +inline double Helix::phase() const {return mPhase;} -inline double Helix::x(double s) const { - if (mSingularity) - return mOrigin.x - s * mCosDipAngle * mSinPhase; - else - return mOrigin.x + (cos(mPhase + s * mH * mCurvature * mCosDipAngle) - mCosPhase) / mCurvature; +inline double Helix::x(double s) const +{ + if (mSingularity) + return mOrigin.x - s*mCosDipAngle*mSinPhase; + else + return mOrigin.x + (cos(mPhase + s*mH*mCurvature*mCosDipAngle)-mCosPhase)/mCurvature; } - -inline double Helix::y(double s) const { - if (mSingularity) - return mOrigin.y + s * mCosDipAngle * mCosPhase; - else - return mOrigin.y + (sin(mPhase + s * mH * mCurvature * mCosDipAngle) - mSinPhase) / mCurvature; + +inline double Helix::y(double s) const +{ + if (mSingularity) + return mOrigin.y + s*mCosDipAngle*mCosPhase; + else + return mOrigin.y + (sin(mPhase + s*mH*mCurvature*mCosDipAngle)-mSinPhase)/mCurvature; } -inline double Helix::z(double s) const { return mOrigin.z + s * mSinDipAngle; } +inline double Helix::z(double s) const +{ + return mOrigin.z + s*mSinDipAngle; +} -inline double Helix::cx(double s) const { - if (mSingularity) - return -mCosDipAngle * mSinPhase; - else - return -sin(mPhase + s * mH * mCurvature * mCosDipAngle) * mH * mCosDipAngle; +inline double Helix::cx(double s) const +{ + if (mSingularity) + return -mCosDipAngle*mSinPhase; + else + return -sin(mPhase + s*mH*mCurvature*mCosDipAngle)*mH*mCosDipAngle; } -inline double Helix::cy(double s) const { - if (mSingularity) - return mCosDipAngle * mCosPhase; - else - return cos(mPhase + s * mH * mCurvature * mCosDipAngle) * mH * mCosDipAngle; +inline double Helix::cy(double s) const +{ + if (mSingularity) + return mCosDipAngle*mCosPhase; + else + return cos(mPhase + s*mH*mCurvature*mCosDipAngle)*mH*mCosDipAngle; } -inline double Helix::cz(double /* s */) const { return mSinDipAngle; } +inline double Helix::cz(double /* s */) const +{ + return mSinDipAngle; +} -inline const edm4hep::Vector3f& Helix::origin() const { return mOrigin; } +inline const edm4hep::Vector3f& Helix::origin() const {return mOrigin;} -inline edm4hep::Vector3f Helix::at(double s) const { return edm4hep::Vector3f(x(s), y(s), z(s)); } +inline edm4hep::Vector3f Helix::at(double s) const +{ + return edm4hep::Vector3f(x(s), y(s), z(s)); +} -inline edm4hep::Vector3f Helix::cat(double s) const { - return edm4hep::Vector3f(cx(s), cy(s), cz(s)); +inline edm4hep::Vector3f Helix::cat(double s) const +{ + return edm4hep::Vector3f(cx(s), cy(s), cz(s)); } -inline double Helix::pathLength(double X, double Y) const { - return fudgePathLength(edm4hep::Vector3f(X, Y, 0)); +inline double Helix::pathLength(double X, double Y) const +{ + return fudgePathLength(edm4hep::Vector3f(X, Y, 0)); } -inline int Helix::bad(double WorldSize) const { +inline int Helix::bad(double WorldSize) const +{ - // int ierr; - if (!::finite(mDipAngle)) - return 11; - if (!::finite(mCurvature)) - return 12; +// int ierr; + if (!::finite(mDipAngle )) return 11; + if (!::finite(mCurvature )) return 12; - // ierr = mOrigin.bad(WorldSize); - // if (ierr) return 3+ierr*100; +// ierr = mOrigin.bad(WorldSize); +// if (ierr) return 3+ierr*100; - if (::fabs(mDipAngle) > 1.58) - return 21; - double qwe = ::fabs(::fabs(mDipAngle) - M_PI / 2); - if (qwe < 1. / WorldSize) - return 31; + if (::fabs(mDipAngle) >1.58) return 21; + double qwe = ::fabs(::fabs(mDipAngle)-M_PI/2); + if (qwe < 1./WorldSize ) return 31; - if (::fabs(mCurvature) > WorldSize) - return 22; - if (mCurvature < 0) - return 32; + if (::fabs(mCurvature) > WorldSize) return 22; + if (mCurvature < 0 ) return 32; - if (abs(mH) != 1) - return 24; + if (abs(mH) != 1 ) return 24; - return 0; + return 0; } -} // namespace eicrecon +} From 4819ae165a88151ee09863b777eb93d583654c9c Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 20 Oct 2025 22:58:23 +0000 Subject: [PATCH 11/17] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/algorithms/reco/Helix.cc | 1029 ++++++++++++++++------------------ src/algorithms/reco/Helix.h | 338 ++++++----- 2 files changed, 658 insertions(+), 709 deletions(-) diff --git a/src/algorithms/reco/Helix.cc b/src/algorithms/reco/Helix.cc index a675330c1f..5534bc31ca 100644 --- a/src/algorithms/reco/Helix.cc +++ b/src/algorithms/reco/Helix.cc @@ -25,15 +25,13 @@ namespace eicrecon { const double Helix::NoSolution = 3.e+33; // basic constructor -Helix::Helix(double c, double d, double phase, const edm4hep::Vector3f& o, int h) -{ - setParameters(c, d, phase, o, h); +Helix::Helix(double c, double d, double phase, const edm4hep::Vector3f& o, int h) { + setParameters(c, d, phase, o, h); } // momentum in GeV, position in cm, B in Tesla -Helix::Helix(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q) -{ - setParameters(p, o, B, q); +Helix::Helix(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q) { + setParameters(p, o, B, q); } // construct using ReconstructParticle @@ -50,508 +48,481 @@ Helix::Helix(const edm4eic::ReconstructedParticle& p, const double b_field) { // construct using TrackParameteters void Helix::setParameters(const edm4eic::TrackParameters& trk, const double b_field) { - const auto mom = edm4hep::utils::sphericalToVector(1.0 / std::abs(trk.getQOverP()), trk.getTheta(), trk.getPhi()); + const auto mom = edm4hep::utils::sphericalToVector(1.0 / std::abs(trk.getQOverP()), + trk.getTheta(), trk.getPhi()); const auto charge = std::copysign(1., trk.getQOverP()); - const auto phi = trk.getPhi(); - const auto loc = trk.getLoc(); - edm4hep::Vector3f pos( -1. * loc.a * sin(phi), loc.a * cos(phi), loc.b); // PCA point + const auto phi = trk.getPhi(); + const auto loc = trk.getLoc(); + edm4hep::Vector3f pos(-1. * loc.a * sin(phi), loc.a * cos(phi), loc.b); // PCA point - setParameters(mom*edm4eic::unit::GeV/dd4hep::GeV, pos*edm4eic::unit::mm/edm4eic::unit::cm, b_field, charge); + setParameters(mom * edm4eic::unit::GeV / dd4hep::GeV, pos * edm4eic::unit::mm / edm4eic::unit::cm, + b_field, charge); } -void Helix::setParameters(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q) -{ - mH = (q*B <= 0) ? 1 : -1; - if(p.y == 0 && p.x == 0) - setPhase((M_PI/4)*(1-2.*mH)); - else - setPhase(atan2(p.y,p.x)-mH*M_PI/2); - setDipAngle(atan2(p.z,edm4hep::utils::magnitudeTransverse(p))); - mOrigin = o; - - double curvature_val = fabs((dd4hep::c_light*dd4hep::nanosecond/dd4hep::meter*q*B/dd4hep::tesla)/ - (edm4hep::utils::magnitude(p)/dd4hep::GeV*mCosDipAngle)/dd4hep::meter); - setCurvature(curvature_val); +void Helix::setParameters(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, + const int q) { + mH = (q * B <= 0) ? 1 : -1; + if (p.y == 0 && p.x == 0) + setPhase((M_PI / 4) * (1 - 2. * mH)); + else + setPhase(atan2(p.y, p.x) - mH * M_PI / 2); + setDipAngle(atan2(p.z, edm4hep::utils::magnitudeTransverse(p))); + mOrigin = o; + + double curvature_val = + fabs((dd4hep::c_light * dd4hep::nanosecond / dd4hep::meter * q * B / dd4hep::tesla) / + (edm4hep::utils::magnitude(p) / dd4hep::GeV * mCosDipAngle) / dd4hep::meter); + setCurvature(curvature_val); } -void Helix::setParameters(double c, double dip, double phase, - const edm4hep::Vector3f& o, int h) -{ - // - // The order in which the parameters are set is important - // since setCurvature might have to adjust the others. - // - mH = (h>=0) ? 1 : -1; // Default is: positive particle - // positive field - mOrigin = o; - setDipAngle(dip); - setPhase(phase); - - // - // Check for singularity and correct for negative curvature. - // May change mH and mPhase. Must therefore be set last. - // - setCurvature(c); - - // - // For the case B=0, h is ill defined. In the following we - // always assume h = +1. Since phase = psi - h * pi/2 - // we have to correct the phase in case h = -1. - // This assumes that the user uses the same h for phase - // as the one he passed to the constructor. - // - if (mSingularity && mH == -1) { - mH = +1; - setPhase(mPhase-M_PI); - } +void Helix::setParameters(double c, double dip, double phase, const edm4hep::Vector3f& o, int h) { + // + // The order in which the parameters are set is important + // since setCurvature might have to adjust the others. + // + mH = (h >= 0) ? 1 : -1; // Default is: positive particle + // positive field + mOrigin = o; + setDipAngle(dip); + setPhase(phase); + + // + // Check for singularity and correct for negative curvature. + // May change mH and mPhase. Must therefore be set last. + // + setCurvature(c); + + // + // For the case B=0, h is ill defined. In the following we + // always assume h = +1. Since phase = psi - h * pi/2 + // we have to correct the phase in case h = -1. + // This assumes that the user uses the same h for phase + // as the one he passed to the constructor. + // + if (mSingularity && mH == -1) { + mH = +1; + setPhase(mPhase - M_PI); + } } -void Helix::setCurvature(double val) -{ - if (val < 0) { - mCurvature = -val; - mH = -mH; - setPhase(mPhase+M_PI); - } - else - mCurvature = val; - - if (fabs(mCurvature) <= std::numeric_limits::epsilon()) - mSingularity = true; // straight line - else - mSingularity = false; // curved +void Helix::setCurvature(double val) { + if (val < 0) { + mCurvature = -val; + mH = -mH; + setPhase(mPhase + M_PI); + } else + mCurvature = val; + + if (fabs(mCurvature) <= std::numeric_limits::epsilon()) + mSingularity = true; // straight line + else + mSingularity = false; // curved } -void Helix::setPhase(double val) -{ - mPhase = val; - mCosPhase = cos(mPhase); - mSinPhase = sin(mPhase); - if (fabs(mPhase) > M_PI) - mPhase = atan2(mSinPhase, mCosPhase); // force range [-pi,pi] +void Helix::setPhase(double val) { + mPhase = val; + mCosPhase = cos(mPhase); + mSinPhase = sin(mPhase); + if (fabs(mPhase) > M_PI) + mPhase = atan2(mSinPhase, mCosPhase); // force range [-pi,pi] } -void Helix::setDipAngle(double val) -{ - mDipAngle = val; - mCosDipAngle = cos(mDipAngle); - mSinDipAngle = sin(mDipAngle); +void Helix::setDipAngle(double val) { + mDipAngle = val; + mCosDipAngle = cos(mDipAngle); + mSinDipAngle = sin(mDipAngle); } -double Helix::xcenter() const -{ - if (mSingularity) - return 0; - else - return mOrigin.x-mCosPhase/mCurvature; +double Helix::xcenter() const { + if (mSingularity) + return 0; + else + return mOrigin.x - mCosPhase / mCurvature; } -double Helix::ycenter() const -{ - if (mSingularity) - return 0; - else - return mOrigin.y-mSinPhase/mCurvature; +double Helix::ycenter() const { + if (mSingularity) + return 0; + else + return mOrigin.y - mSinPhase / mCurvature; } -double Helix::fudgePathLength(const edm4hep::Vector3f& p) const -{ - double s; - double dx = p.x-mOrigin.x; - double dy = p.y-mOrigin.y; - - if (mSingularity) { - s = (dy*mCosPhase - dx*mSinPhase)/mCosDipAngle; - } - else { - s = atan2(dy*mCosPhase - dx*mSinPhase, - 1/mCurvature + dx*mCosPhase+dy*mSinPhase)/ - (mH*mCurvature*mCosDipAngle); - } - return s; +double Helix::fudgePathLength(const edm4hep::Vector3f& p) const { + double s; + double dx = p.x - mOrigin.x; + double dy = p.y - mOrigin.y; + + if (mSingularity) { + s = (dy * mCosPhase - dx * mSinPhase) / mCosDipAngle; + } else { + s = atan2(dy * mCosPhase - dx * mSinPhase, 1 / mCurvature + dx * mCosPhase + dy * mSinPhase) / + (mH * mCurvature * mCosDipAngle); + } + return s; } -double Helix::distance(const edm4hep::Vector3f& p, bool scanPeriods) const -{ - return edm4hep::utils::magnitude(this->at(pathLength(p,scanPeriods))-p); +double Helix::distance(const edm4hep::Vector3f& p, bool scanPeriods) const { + return edm4hep::utils::magnitude(this->at(pathLength(p, scanPeriods)) - p); } -double Helix::pathLength(const edm4hep::Vector3f& p, bool scanPeriods) const -{ +double Helix::pathLength(const edm4hep::Vector3f& p, bool scanPeriods) const { + // + // Returns the path length at the distance of closest + // approach between the helix and point p. + // For the case of B=0 (straight line) the path length + // can be calculated analytically. For B>0 there is + // unfortunately no easy solution to the problem. + // Here we use the Newton method to find the root of the + // referring equation. The 'fudgePathLength' serves + // as a starting value. + // + double s; + double dx = p.x - mOrigin.x; + double dy = p.y - mOrigin.y; + double dz = p.z - mOrigin.z; + + if (mSingularity) { + s = mCosDipAngle * (mCosPhase * dy - mSinPhase * dx) + mSinDipAngle * dz; + } else { // + const double MaxPrecisionNeeded = edm4eic::unit::um; + const int MaxIterations = 100; + // - // Returns the path length at the distance of closest - // approach between the helix and point p. - // For the case of B=0 (straight line) the path length - // can be calculated analytically. For B>0 there is - // unfortunately no easy solution to the problem. - // Here we use the Newton method to find the root of the - // referring equation. The 'fudgePathLength' serves - // as a starting value. + // The math is taken from Maple with C(expr,optimized) and + // some hand-editing. It is not very nice but efficient. // - double s; - double dx = p.x-mOrigin.x; - double dy = p.y-mOrigin.y; - double dz = p.z-mOrigin.z; + double t34 = mCurvature * mCosDipAngle * mCosDipAngle; + double t41 = mSinDipAngle * mSinDipAngle; + double t6, t7, t11, t12, t19; - if (mSingularity) { - s = mCosDipAngle*(mCosPhase*dy-mSinPhase*dx) + - mSinDipAngle*dz; + // + // Get a first guess by using the dca in 2D. Since + // in some extreme cases we might be off by n periods + // we add (subtract) periods in case we get any closer. + // + s = fudgePathLength(p); + + if (scanPeriods) { + double ds = period(); + int j, jmin = 0; + double d, dmin = edm4hep::utils::magnitude(at(s) - p); + for (j = 1; j < MaxIterations; j++) { + if ((d = edm4hep::utils::magnitude(at(s + j * ds) - p)) < dmin) { + dmin = d; + jmin = j; + } else + break; + } + for (j = -1; -j < MaxIterations; j--) { + if ((d = edm4hep::utils::magnitude(at(s + j * ds) - p)) < dmin) { + dmin = d; + jmin = j; + } else + break; + } + if (jmin) + s += jmin * ds; } - else { // - const double MaxPrecisionNeeded = edm4eic::unit::um; - const int MaxIterations = 100; - - // - // The math is taken from Maple with C(expr,optimized) and - // some hand-editing. It is not very nice but efficient. - // - double t34 = mCurvature*mCosDipAngle*mCosDipAngle; - double t41 = mSinDipAngle*mSinDipAngle; - double t6, t7, t11, t12, t19; - - // - // Get a first guess by using the dca in 2D. Since - // in some extreme cases we might be off by n periods - // we add (subtract) periods in case we get any closer. - // - s = fudgePathLength(p); - - if (scanPeriods) { - double ds = period(); - int j, jmin = 0; - double d, dmin = edm4hep::utils::magnitude(at(s) - p); - for(j=1; j::max(); - else - return fabs(2*M_PI/(mH*mCurvature*mCosDipAngle)); +double Helix::period() const { + if (mSingularity) + return std::numeric_limits::max(); + else + return fabs(2 * M_PI / (mH * mCurvature * mCosDipAngle)); } -std::pair Helix::pathLength(double r) const -{ - std::pair value; - std::pair VALUE(999999999.,999999999.); +std::pair Helix::pathLength(double r) const { + std::pair value; + std::pair VALUE(999999999., 999999999.); + // + // The math is taken from Maple with C(expr,optimized) and + // some hand-editing. It is not very nice but efficient. + // 'first' is the smallest of the two solutions (may be negative) + // 'second' is the other. + // + if (mSingularity) { + double t1 = mCosDipAngle * (mOrigin.x * mSinPhase - mOrigin.y * mCosPhase); + double t12 = mOrigin.y * mOrigin.y; + double t13 = mCosPhase * mCosPhase; + double t15 = r * r; + double t16 = mOrigin.x * mOrigin.x; + double t20 = + -mCosDipAngle * mCosDipAngle * + (2.0 * mOrigin.x * mSinPhase * mOrigin.y * mCosPhase + t12 - t12 * t13 - t15 + t13 * t16); + if (t20 < 0.) + return VALUE; + t20 = ::sqrt(t20); + value.first = (t1 - t20) / (mCosDipAngle * mCosDipAngle); + value.second = (t1 + t20) / (mCosDipAngle * mCosDipAngle); + } else { + double t1 = mOrigin.y * mCurvature; + double t2 = mSinPhase; + double t3 = mCurvature * mCurvature; + double t4 = mOrigin.y * t2; + double t5 = mCosPhase; + double t6 = mOrigin.x * t5; + double t8 = mOrigin.x * mOrigin.x; + double t11 = mOrigin.y * mOrigin.y; + double t14 = r * r; + double t15 = t14 * mCurvature; + double t17 = t8 * t8; + double t19 = t11 * t11; + double t21 = t11 * t3; + double t23 = t5 * t5; + double t32 = t14 * t14; + double t35 = t14 * t3; + double t38 = 8.0 * t4 * t6 - 4.0 * t1 * t2 * t8 - 4.0 * t11 * mCurvature * t6 + 4.0 * t15 * t6 + + t17 * t3 + t19 * t3 + 2.0 * t21 * t8 + 4.0 * t8 * t23 - + 4.0 * t8 * mOrigin.x * mCurvature * t5 - 4.0 * t11 * t23 - + 4.0 * t11 * mOrigin.y * mCurvature * t2 + 4.0 * t11 - 4.0 * t14 + t32 * t3 + + 4.0 * t15 * t4 - 2.0 * t35 * t11 - 2.0 * t35 * t8; + double t40 = (-t3 * t38); + if (t40 < 0.) + return VALUE; + t40 = ::sqrt(t40); + + double t43 = mOrigin.x * mCurvature; + double t45 = 2.0 * t5 - t35 + t21 + 2.0 - 2.0 * t1 * t2 - 2.0 * t43 - 2.0 * t43 * t5 + t8 * t3; + double t46 = mH * mCosDipAngle * mCurvature; + + value.first = (-mPhase + 2.0 * atan((-2.0 * t1 + 2.0 * t2 + t40) / t45)) / t46; + value.second = -(mPhase + 2.0 * atan((2.0 * t1 - 2.0 * t2 + t40) / t45)) / t46; + // - // The math is taken from Maple with C(expr,optimized) and - // some hand-editing. It is not very nice but efficient. - // 'first' is the smallest of the two solutions (may be negative) - // 'second' is the other. + // Solution can be off by +/- one period, select smallest // - if (mSingularity) { - double t1 = mCosDipAngle*(mOrigin.x*mSinPhase-mOrigin.y*mCosPhase); - double t12 = mOrigin.y*mOrigin.y; - double t13 = mCosPhase*mCosPhase; - double t15 = r*r; - double t16 = mOrigin.x*mOrigin.x; - double t20 = -mCosDipAngle*mCosDipAngle*(2.0*mOrigin.x*mSinPhase*mOrigin.y*mCosPhase + - t12-t12*t13-t15+t13*t16); - if (t20<0.) return VALUE; - t20 = ::sqrt(t20); - value.first = (t1-t20)/(mCosDipAngle*mCosDipAngle); - value.second = (t1+t20)/(mCosDipAngle*mCosDipAngle); + double p = period(); + if (!std::isnan(value.first)) { + if (fabs(value.first - p) < fabs(value.first)) + value.first = value.first - p; + else if (fabs(value.first + p) < fabs(value.first)) + value.first = value.first + p; } - else { - double t1 = mOrigin.y*mCurvature; - double t2 = mSinPhase; - double t3 = mCurvature*mCurvature; - double t4 = mOrigin.y*t2; - double t5 = mCosPhase; - double t6 = mOrigin.x*t5; - double t8 = mOrigin.x*mOrigin.x; - double t11 = mOrigin.y*mOrigin.y; - double t14 = r*r; - double t15 = t14*mCurvature; - double t17 = t8*t8; - double t19 = t11*t11; - double t21 = t11*t3; - double t23 = t5*t5; - double t32 = t14*t14; - double t35 = t14*t3; - double t38 = 8.0*t4*t6 - 4.0*t1*t2*t8 - 4.0*t11*mCurvature*t6 + - 4.0*t15*t6 + t17*t3 + t19*t3 + 2.0*t21*t8 + 4.0*t8*t23 - - 4.0*t8*mOrigin.x*mCurvature*t5 - 4.0*t11*t23 - - 4.0*t11*mOrigin.y*mCurvature*t2 + 4.0*t11 - 4.0*t14 + - t32*t3 + 4.0*t15*t4 - 2.0*t35*t11 - 2.0*t35*t8; - double t40 = (-t3*t38); - if (t40<0.) return VALUE; - t40 = ::sqrt(t40); - - double t43 = mOrigin.x*mCurvature; - double t45 = 2.0*t5 - t35 + t21 + 2.0 - 2.0*t1*t2 -2.0*t43 - 2.0*t43*t5 + t8*t3; - double t46 = mH*mCosDipAngle*mCurvature; - - value.first = (-mPhase + 2.0*atan((-2.0*t1 + 2.0*t2 + t40)/t45))/t46; - value.second = -(mPhase + 2.0*atan((2.0*t1 - 2.0*t2 + t40)/t45))/t46; - - // - // Solution can be off by +/- one period, select smallest - // - double p = period(); - if (! std::isnan(value.first)) { - if (fabs(value.first-p) < fabs(value.first)) value.first = value.first-p; - else if (fabs(value.first+p) < fabs(value.first)) value.first = value.first+p; - } - if (! std::isnan(value.second)) { - if (fabs(value.second-p) < fabs(value.second)) value.second = value.second-p; - else if (fabs(value.second+p) < fabs(value.second)) value.second = value.second+p; - } + if (!std::isnan(value.second)) { + if (fabs(value.second - p) < fabs(value.second)) + value.second = value.second - p; + else if (fabs(value.second + p) < fabs(value.second)) + value.second = value.second + p; } - if (value.first > value.second) - std::swap(value.first,value.second); - return(value); + } + if (value.first > value.second) + std::swap(value.first, value.second); + return (value); } -std::pair Helix::pathLength(double r, double x, double y) -{ - double x0 = mOrigin.x; - double y0 = mOrigin.y; - mOrigin.x = x0-x; - mOrigin.y = y0-y; - std::pair result = this->pathLength(r); - mOrigin.x = x0; - mOrigin.y = y0; - return result; +std::pair Helix::pathLength(double r, double x, double y) { + double x0 = mOrigin.x; + double y0 = mOrigin.y; + mOrigin.x = x0 - x; + mOrigin.y = y0 - y; + std::pair result = this->pathLength(r); + mOrigin.x = x0; + mOrigin.y = y0; + return result; } -double Helix::pathLength(const edm4hep::Vector3f& r, - const edm4hep::Vector3f& n) const -{ +double Helix::pathLength(const edm4hep::Vector3f& r, const edm4hep::Vector3f& n) const { + // + // Vector 'r' defines the position of the center and + // vector 'n' the normal vector of the plane. + // For a straight line there is a simple analytical + // solution. For curvatures > 0 the root is determined + // by Newton method. In case no valid s can be found + // the max. largest value for s is returned. + // + double s; + + if (mSingularity) { + double t = n.z * mSinDipAngle + n.y * mCosDipAngle * mCosPhase - n.x * mCosDipAngle * mSinPhase; + if (t == 0) + s = NoSolution; + else + s = ((r - mOrigin) * n) / t; + } else { + const double MaxPrecisionNeeded = edm4eic::unit::um; + const int MaxIterations = 20; + + double A = mCurvature * ((mOrigin - r) * n) - n.x * mCosPhase - n.y * mSinPhase; + double t = mH * mCurvature * mCosDipAngle; + double u = n.z * mCurvature * mSinDipAngle; + + double a, f, fp; + double sOld = s = 0; + double shiftOld = 0; + double shift; + // (cos(angMax)-1)/angMax = 0.1 + const double angMax = 0.21; + double deltas = fabs(angMax / (mCurvature * mCosDipAngle)); + // dampingFactor = exp(-0.5); + // double dampingFactor = 0.60653; + int i; + + for (i = 0; i < MaxIterations; i++) { + a = t * s + mPhase; + double sina = sin(a); + double cosa = cos(a); + f = A + n.x * cosa + n.y * sina + u * s; + fp = -n.x * sina * t + n.y * cosa * t + u; + if (fabs(fp) * deltas <= fabs(f)) { //too big step + int sgn = 1; + if (fp < 0.) + sgn = -sgn; + if (f < 0.) + sgn = -sgn; + shift = sgn * deltas; + if (shift < 0) + shift *= 0.9; // don't get stuck shifting +/-deltas + } else { + shift = f / fp; + } + s -= shift; + shiftOld = shift; + if (fabs(sOld - s) < MaxPrecisionNeeded) + break; + sOld = s; + } + if (i == MaxIterations) + return NoSolution; + } + return s; +} + +std::pair Helix::pathLengths(const Helix& h, double minStepSize, + double minRange) const { + // + // Cannot handle case where one is a helix + // and the other one is a straight line. + // + if (mSingularity != h.mSingularity) + return std::pair(NoSolution, NoSolution); + + double s1, s2; + + if (mSingularity) { // - // Vector 'r' defines the position of the center and - // vector 'n' the normal vector of the plane. - // For a straight line there is a simple analytical - // solution. For curvatures > 0 the root is determined - // by Newton method. In case no valid s can be found - // the max. largest value for s is returned. + // Analytic solution // - double s; + edm4hep::Vector3f dv = h.mOrigin - mOrigin; + edm4hep::Vector3f a(-mCosDipAngle * mSinPhase, mCosDipAngle * mCosPhase, mSinDipAngle); + edm4hep::Vector3f b(-h.mCosDipAngle * h.mSinPhase, h.mCosDipAngle * h.mCosPhase, + h.mSinDipAngle); + double ab = a * b; + double g = dv * a; + double k = dv * b; + s2 = (k - ab * g) / (ab * ab - 1.); + s1 = g + s2 * ab; + return std::pair(s1, s2); + } else { + // + // First step: get dca in the xy-plane as start value + // + double dx = h.xcenter() - xcenter(); + double dy = h.ycenter() - ycenter(); + double dd = ::sqrt(dx * dx + dy * dy); + double r1 = 1 / curvature(); + double r2 = 1 / h.curvature(); - if (mSingularity) { - double t = n.z*mSinDipAngle + - n.y*mCosDipAngle*mCosPhase - - n.x*mCosDipAngle*mSinPhase; - if (t == 0) - s = NoSolution; - else - s = ((r - mOrigin)*n)/t; - } - else { - const double MaxPrecisionNeeded = edm4eic::unit::um; - const int MaxIterations = 20; - - double A = mCurvature*((mOrigin - r)*n) - - n.x*mCosPhase - - n.y*mSinPhase; - double t = mH*mCurvature*mCosDipAngle; - double u = n.z*mCurvature*mSinDipAngle; - - double a, f, fp; - double sOld = s = 0; - double shiftOld = 0; - double shift; -// (cos(angMax)-1)/angMax = 0.1 - const double angMax = 0.21; - double deltas = fabs(angMax/(mCurvature*mCosDipAngle)); -// dampingFactor = exp(-0.5); -// double dampingFactor = 0.60653; - int i; - - for (i=0; i dd ? -1 : 1); // set -1 when *this* helix is + // completely contained in the other + x = xcenter() + rsign * r1 * dx / dd; + y = ycenter() + rsign * r1 * dy / dd; + s = pathLength(x, y); } - return s; -} -std::pair -Helix::pathLengths(const Helix& h, double minStepSize, double minRange) const -{ // - // Cannot handle case where one is a helix - // and the other one is a straight line. + // Second step: scan in decreasing intervals around seed 's' + // minRange and minStepSize are passed as arguments to the method. + // They have default values defined in the header file. // - if (mSingularity != h.mSingularity) - return std::pair(NoSolution, NoSolution); - - double s1, s2; - - if (mSingularity) { - // - // Analytic solution - // - edm4hep::Vector3f dv = h.mOrigin - mOrigin; - edm4hep::Vector3f a(-mCosDipAngle*mSinPhase, - mCosDipAngle*mCosPhase, - mSinDipAngle); - edm4hep::Vector3f b(-h.mCosDipAngle*h.mSinPhase, - h.mCosDipAngle*h.mCosPhase, - h.mSinDipAngle); - double ab = a*b; - double g = dv*a; - double k = dv*b; - s2 = (k-ab*g)/(ab*ab-1.); - s1 = g+s2*ab; - return std::pair(s1, s2); - } - else { - // - // First step: get dca in the xy-plane as start value - // - double dx = h.xcenter() - xcenter(); - double dy = h.ycenter() - ycenter(); - double dd = ::sqrt(dx*dx + dy*dy); - double r1 = 1/curvature(); - double r2 = 1/h.curvature(); - - double cosAlpha = (r1*r1 + dd*dd - r2*r2)/(2*r1*dd); - - double s; - double x, y; - if (fabs(cosAlpha) < 1) { // two solutions - double sinAlpha = sin(acos(cosAlpha)); - x = xcenter() + r1*(cosAlpha*dx - sinAlpha*dy)/dd; - y = ycenter() + r1*(sinAlpha*dx + cosAlpha*dy)/dd; - s = pathLength(x, y); - x = xcenter() + r1*(cosAlpha*dx + sinAlpha*dy)/dd; - y = ycenter() + r1*(cosAlpha*dy - sinAlpha*dx)/dd; - double a = pathLength(x, y); - if (h.distance(at(a)) < h.distance(at(s))) s = a; - } - else { // no intersection (or exactly one) - int rsign = ((r2-r1) > dd ? -1 : 1); // set -1 when *this* helix is - // completely contained in the other - x = xcenter() + rsign*r1*dx/dd; - y = ycenter() + rsign*r1*dy/dd; - s = pathLength(x, y); - } - - // - // Second step: scan in decreasing intervals around seed 's' - // minRange and minStepSize are passed as arguments to the method. - // They have default values defined in the header file. - // - double dmin = h.distance(at(s)); - double range = std::max(2*dmin, minRange); - double ds = range/10; - double slast=-999999, ss, d; - s1 = s - range/2.; - s2 = s + range/2.; - - while (ds > minStepSize) { - for (ss=s1; ss minStepSize) { + for (ss = s1; ss < s2 + ds; ss += ds) { + d = h.distance(at(ss)); + if (d < dmin) { + dmin = d; + s = ss; } - return std::pair(s, h.pathLength(at(s))); + slast = ss; + } + // + // In the rare cases where the minimum is at the + // the border of the current range we shift the range + // and start all over, i.e we do not decrease 'ds'. + // Else we decrease the search interval around the + // current minimum and redo the scan in smaller steps. + // + if (s == s1) { + d = 0.8 * (s2 - s1); + s1 -= d; + s2 -= d; + } else if (s == slast) { + d = 0.8 * (s2 - s1); + s1 += d; + s2 += d; + } else { + s1 = s - ds; + s2 = s + ds; + ds /= 10; + } } + return std::pair(s, h.pathLength(at(s))); + } } - -void Helix::moveOrigin(double s) -{ - if (mSingularity) - mOrigin = at(s); - else { - edm4hep::Vector3f newOrigin = at(s); - double newPhase = atan2(newOrigin.y - ycenter(), - newOrigin.x - xcenter()); - mOrigin = newOrigin; - setPhase(newPhase); - } +void Helix::moveOrigin(double s) { + if (mSingularity) + mOrigin = at(s); + else { + edm4hep::Vector3f newOrigin = at(s); + double newPhase = atan2(newOrigin.y - ycenter(), newOrigin.x - xcenter()); + mOrigin = newOrigin; + setPhase(newPhase); + } } /* int operator== (const Helix& a, const Helix& b) @@ -560,100 +531,88 @@ int operator== (const Helix& a, const Helix& b) // Checks for numerical identity only ! // return (a.origin() == b.origin() && - a.dipAngle() == b.dipAngle() && - a.curvature() == b.curvature() && - a.phase() == b.phase() && - a.h() == b.h()); + a.dipAngle() == b.dipAngle() && + a.curvature() == b.curvature() && + a.phase() == b.phase() && + a.h() == b.h()); } int operator!= (const Helix& a, const Helix& b) {return !(a == b);} */ //std::ostream& operator<<(std::ostream& os, const Helix& h) -void Helix::Print() const -{ - std::cout << "(" - << "curvature = " << mCurvature << ", " - << "dip angle = " << mDipAngle << ", " - << "phase = " << mPhase << ", " - << "h = " << mH << ", " - << "origin = " << mOrigin.x << " " << mOrigin.y << " " << mOrigin.z << ")" - << std::endl; +void Helix::Print() const { + std::cout << "(" + << "curvature = " << mCurvature << ", " + << "dip angle = " << mDipAngle << ", " + << "phase = " << mPhase << ", " + << "h = " << mH << ", " + << "origin = " << mOrigin.x << " " << mOrigin.y << " " << mOrigin.z << ")" << std::endl; } -edm4hep::Vector3f Helix::momentum(double B) const -{ - if (mSingularity) - return(edm4hep::Vector3f(0,0,0)); - else { - double pt = edm4eic::unit::GeV*fabs(dd4hep::c_light*dd4hep::nanosecond/dd4hep::meter*B/dd4hep::tesla)/(fabs(mCurvature)*dd4hep::meter); - - return (edm4hep::Vector3f(pt*cos(mPhase+mH*M_PI/2), // pos part pos field - pt*sin(mPhase+mH*M_PI/2), - pt*tan(mDipAngle))); - } -} +edm4hep::Vector3f Helix::momentum(double B) const { + if (mSingularity) + return (edm4hep::Vector3f(0, 0, 0)); + else { + double pt = edm4eic::unit::GeV * + fabs(dd4hep::c_light * dd4hep::nanosecond / dd4hep::meter * B / dd4hep::tesla) / + (fabs(mCurvature) * dd4hep::meter); -edm4hep::Vector3f Helix::momentumAt(double S, double B) const -{ - // Obtain phase-shifted momentum from phase-shift of origin - Helix tmp(*this); - tmp.moveOrigin(S); - return tmp.momentum(B); + return (edm4hep::Vector3f(pt * cos(mPhase + mH * M_PI / 2), // pos part pos field + pt * sin(mPhase + mH * M_PI / 2), pt * tan(mDipAngle))); + } } -int Helix::charge(double B) const -{ - return (B > 0 ? -mH : mH); +edm4hep::Vector3f Helix::momentumAt(double S, double B) const { + // Obtain phase-shifted momentum from phase-shift of origin + Helix tmp(*this); + tmp.moveOrigin(S); + return tmp.momentum(B); } -double Helix::geometricSignedDistance(double x, double y) -{ - // Geometric signed distance - double thePath = this->pathLength(x,y); - edm4hep::Vector3f DCA2dPosition = this->at(thePath); - DCA2dPosition.z = 0; - edm4hep::Vector3f position(x,y,0); - edm4hep::Vector3f DCAVec = (DCA2dPosition-position); - edm4hep::Vector3f momVec; - // Deal with straight tracks - if (this->mSingularity) { - momVec = this->at(1)- this->at(0); - momVec.z = 0; - } - else { - momVec = this->momentumAt(thePath,1./dd4hep::tesla); // Don't care about Bmag. Helicity is what matters. - momVec.z = 0; - } - - double cross = DCAVec.x*momVec.y - DCAVec.y*momVec.x; - double theSign = (cross>=0) ? 1. : -1.; - return theSign*edm4hep::utils::magnitudeTransverse(DCAVec); +int Helix::charge(double B) const { return (B > 0 ? -mH : mH); } + +double Helix::geometricSignedDistance(double x, double y) { + // Geometric signed distance + double thePath = this->pathLength(x, y); + edm4hep::Vector3f DCA2dPosition = this->at(thePath); + DCA2dPosition.z = 0; + edm4hep::Vector3f position(x, y, 0); + edm4hep::Vector3f DCAVec = (DCA2dPosition - position); + edm4hep::Vector3f momVec; + // Deal with straight tracks + if (this->mSingularity) { + momVec = this->at(1) - this->at(0); + momVec.z = 0; + } else { + momVec = this->momentumAt( + thePath, 1. / dd4hep::tesla); // Don't care about Bmag. Helicity is what matters. + momVec.z = 0; + } + + double cross = DCAVec.x * momVec.y - DCAVec.y * momVec.x; + double theSign = (cross >= 0) ? 1. : -1.; + return theSign * edm4hep::utils::magnitudeTransverse(DCAVec); } -double Helix::curvatureSignedDistance(double x, double y) -{ - // Protect against mH = 0 or zero field - if (this->mSingularity || abs(this->mH)<=0) { - return (this->geometricSignedDistance(x,y)); - } - else { - return (this->geometricSignedDistance(x,y))/(this->mH); - } - +double Helix::curvatureSignedDistance(double x, double y) { + // Protect against mH = 0 or zero field + if (this->mSingularity || abs(this->mH) <= 0) { + return (this->geometricSignedDistance(x, y)); + } else { + return (this->geometricSignedDistance(x, y)) / (this->mH); + } } -double Helix::geometricSignedDistance(const edm4hep::Vector3f& pos) -{ - double sdca2d = this->geometricSignedDistance(pos.x,pos.y); - double theSign = (sdca2d>=0) ? 1. : -1.; - return (this->distance(pos))*theSign; +double Helix::geometricSignedDistance(const edm4hep::Vector3f& pos) { + double sdca2d = this->geometricSignedDistance(pos.x, pos.y); + double theSign = (sdca2d >= 0) ? 1. : -1.; + return (this->distance(pos)) * theSign; } -double Helix::curvatureSignedDistance(const edm4hep::Vector3f& pos) -{ - double sdca2d = this->curvatureSignedDistance(pos.x,pos.y); - double theSign = (sdca2d>=0) ? 1. : -1.; - return (this->distance(pos))*theSign; +double Helix::curvatureSignedDistance(const edm4hep::Vector3f& pos) { + double sdca2d = this->curvatureSignedDistance(pos.x, pos.y); + double theSign = (sdca2d >= 0) ? 1. : -1.; + return (this->distance(pos)) * theSign; } -} // namespace eicrecon \ No newline at end of file +} // namespace eicrecon diff --git a/src/algorithms/reco/Helix.h b/src/algorithms/reco/Helix.h index 437336343d..8e7f710c60 100644 --- a/src/algorithms/reco/Helix.h +++ b/src/algorithms/reco/Helix.h @@ -21,119 +21,118 @@ namespace eicrecon { class Helix { - bool mSingularity; // true for straight line case (B=0) - edm4hep::Vector3f mOrigin; - double mDipAngle; - double mCurvature; - double mPhase; - int mH; // -sign(q*B); - - double mCosDipAngle; - double mSinDipAngle; - double mCosPhase; - double mSinPhase; + bool mSingularity; // true for straight line case (B=0) + edm4hep::Vector3f mOrigin; + double mDipAngle; + double mCurvature; + double mPhase; + int mH; // -sign(q*B); + + double mCosDipAngle; + double mSinDipAngle; + double mCosPhase; + double mSinPhase; public: - /// curvature, dip angle, phase, origin, h - Helix(const double c, const double dip, const double phase, const edm4hep::Vector3f& o, const int h=-1); - - /// momentum, origin, b_field, charge - Helix(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q); - - /// ReconstructParticle, b field - Helix(const edm4eic::ReconstructedParticle& p, const double b_field); - - ~Helix() = default; - - double dipAngle() const; - double curvature() const; /// 1/R in xy-plane - double phase() const; /// aziumth in xy-plane measured from ring center - double xcenter() const; /// x-center of circle in xy-plane - double ycenter() const; /// y-center of circle in xy-plane - int h() const; /// -sign(q*B); - - const edm4hep::Vector3f& origin() const; /// starting point - - void setParameters(double c, double dip, double phase, const edm4hep::Vector3f& o, int h); - - void setParameters(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q); - - /// edm4eic::TrackParameters, b field - void setParameters(const edm4eic::TrackParameters& trk, const double b_field); - - /// coordinates of helix at point s - double x(double s) const; - double y(double s) const; - double z(double s) const; - - edm4hep::Vector3f at(double s) const; - - /// pointing vector of helix at point s - double cx(double s) const; - double cy(double s) const; - double cz(double s = 0) const; - - edm4hep::Vector3f cat(double s) const; - - /// returns period length of helix - double period() const; - - /// path length at given r (cylindrical r) - std::pair pathLength(double r) const; - - /// path length at given r (cylindrical r, cylinder axis at x,y) - std::pair pathLength(double r, double x, double y); - - /// path length at distance of closest approach to a given point - double pathLength(const edm4hep::Vector3f& p, bool scanPeriods = true) const; - - /// path length at intersection with plane - double pathLength(const edm4hep::Vector3f& r, - const edm4hep::Vector3f& n) const; - - /// path length at distance of closest approach in the xy-plane to a given point - double pathLength(double x, double y) const; - - /// path lengths at dca between two helices - std::pair pathLengths(const Helix&, - double minStepSize = 10*edm4eic::unit::um, - double minRange = 10*edm4eic::unit::cm) const; - - /// minimal distance between point and helix - double distance(const edm4hep::Vector3f& p, bool scanPeriods = true) const; - - /// checks for valid parametrization - bool valid(double world = 1.e+5) const {return !bad(world);} - int bad(double world = 1.e+5) const; - - /// move the origin along the helix to s which becomes then s=0 - virtual void moveOrigin(double s); - - static const double NoSolution; - - - void setCurvature(double); /// performs also various checks - void setPhase(double); - void setDipAngle(double); - - /// value of S where distance in x-y plane is minimal - double fudgePathLength(const edm4hep::Vector3f&) const; - - // Requires: signed Magnetic Field - edm4hep::Vector3f momentum(double) const; // returns the momentum at origin - edm4hep::Vector3f momentumAt(double, double) const; // returns momentum at S - int charge(double) const; // returns charge of particle - // 2d DCA to x,y point signed relative to curvature - double curvatureSignedDistance(double x, double y) ; - // 2d DCA to x,y point signed relative to rotation - double geometricSignedDistance(double x, double y) ; - // 3d DCA to 3d point signed relative to curvature - double curvatureSignedDistance(const edm4hep::Vector3f&) ; - // 3d DCA to 3d point signed relative to rotation - double geometricSignedDistance(const edm4hep::Vector3f&) ; - - // - void Print() const; + /// curvature, dip angle, phase, origin, h + Helix(const double c, const double dip, const double phase, const edm4hep::Vector3f& o, + const int h = -1); + + /// momentum, origin, b_field, charge + Helix(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q); + + /// ReconstructParticle, b field + Helix(const edm4eic::ReconstructedParticle& p, const double b_field); + + ~Helix() = default; + + double dipAngle() const; + double curvature() const; /// 1/R in xy-plane + double phase() const; /// aziumth in xy-plane measured from ring center + double xcenter() const; /// x-center of circle in xy-plane + double ycenter() const; /// y-center of circle in xy-plane + int h() const; /// -sign(q*B); + + const edm4hep::Vector3f& origin() const; /// starting point + + void setParameters(double c, double dip, double phase, const edm4hep::Vector3f& o, int h); + + void setParameters(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, + const int q); + + /// edm4eic::TrackParameters, b field + void setParameters(const edm4eic::TrackParameters& trk, const double b_field); + + /// coordinates of helix at point s + double x(double s) const; + double y(double s) const; + double z(double s) const; + + edm4hep::Vector3f at(double s) const; + + /// pointing vector of helix at point s + double cx(double s) const; + double cy(double s) const; + double cz(double s = 0) const; + + edm4hep::Vector3f cat(double s) const; + + /// returns period length of helix + double period() const; + + /// path length at given r (cylindrical r) + std::pair pathLength(double r) const; + + /// path length at given r (cylindrical r, cylinder axis at x,y) + std::pair pathLength(double r, double x, double y); + + /// path length at distance of closest approach to a given point + double pathLength(const edm4hep::Vector3f& p, bool scanPeriods = true) const; + + /// path length at intersection with plane + double pathLength(const edm4hep::Vector3f& r, const edm4hep::Vector3f& n) const; + + /// path length at distance of closest approach in the xy-plane to a given point + double pathLength(double x, double y) const; + + /// path lengths at dca between two helices + std::pair pathLengths(const Helix&, double minStepSize = 10 * edm4eic::unit::um, + double minRange = 10 * edm4eic::unit::cm) const; + + /// minimal distance between point and helix + double distance(const edm4hep::Vector3f& p, bool scanPeriods = true) const; + + /// checks for valid parametrization + bool valid(double world = 1.e+5) const { return !bad(world); } + int bad(double world = 1.e+5) const; + + /// move the origin along the helix to s which becomes then s=0 + virtual void moveOrigin(double s); + + static const double NoSolution; + + void setCurvature(double); /// performs also various checks + void setPhase(double); + void setDipAngle(double); + + /// value of S where distance in x-y plane is minimal + double fudgePathLength(const edm4hep::Vector3f&) const; + + // Requires: signed Magnetic Field + edm4hep::Vector3f momentum(double) const; // returns the momentum at origin + edm4hep::Vector3f momentumAt(double, double) const; // returns momentum at S + int charge(double) const; // returns charge of particle + // 2d DCA to x,y point signed relative to curvature + double curvatureSignedDistance(double x, double y); + // 2d DCA to x,y point signed relative to rotation + double geometricSignedDistance(double x, double y); + // 3d DCA to 3d point signed relative to curvature + double curvatureSignedDistance(const edm4hep::Vector3f&); + // 3d DCA to 3d point signed relative to rotation + double geometricSignedDistance(const edm4hep::Vector3f&); + + // + void Print() const; }; // end class Helix @@ -147,92 +146,83 @@ class Helix { // // Inline functions // -inline int Helix::h() const {return mH;} +inline int Helix::h() const { return mH; } -inline double Helix::dipAngle() const {return mDipAngle;} +inline double Helix::dipAngle() const { return mDipAngle; } -inline double Helix::curvature() const {return mCurvature;} +inline double Helix::curvature() const { return mCurvature; } -inline double Helix::phase() const {return mPhase;} +inline double Helix::phase() const { return mPhase; } -inline double Helix::x(double s) const -{ - if (mSingularity) - return mOrigin.x - s*mCosDipAngle*mSinPhase; - else - return mOrigin.x + (cos(mPhase + s*mH*mCurvature*mCosDipAngle)-mCosPhase)/mCurvature; -} - -inline double Helix::y(double s) const -{ - if (mSingularity) - return mOrigin.y + s*mCosDipAngle*mCosPhase; - else - return mOrigin.y + (sin(mPhase + s*mH*mCurvature*mCosDipAngle)-mSinPhase)/mCurvature; +inline double Helix::x(double s) const { + if (mSingularity) + return mOrigin.x - s * mCosDipAngle * mSinPhase; + else + return mOrigin.x + (cos(mPhase + s * mH * mCurvature * mCosDipAngle) - mCosPhase) / mCurvature; } -inline double Helix::z(double s) const -{ - return mOrigin.z + s*mSinDipAngle; +inline double Helix::y(double s) const { + if (mSingularity) + return mOrigin.y + s * mCosDipAngle * mCosPhase; + else + return mOrigin.y + (sin(mPhase + s * mH * mCurvature * mCosDipAngle) - mSinPhase) / mCurvature; } -inline double Helix::cx(double s) const -{ - if (mSingularity) - return -mCosDipAngle*mSinPhase; - else - return -sin(mPhase + s*mH*mCurvature*mCosDipAngle)*mH*mCosDipAngle; +inline double Helix::z(double s) const { return mOrigin.z + s * mSinDipAngle; } + +inline double Helix::cx(double s) const { + if (mSingularity) + return -mCosDipAngle * mSinPhase; + else + return -sin(mPhase + s * mH * mCurvature * mCosDipAngle) * mH * mCosDipAngle; } -inline double Helix::cy(double s) const -{ - if (mSingularity) - return mCosDipAngle*mCosPhase; - else - return cos(mPhase + s*mH*mCurvature*mCosDipAngle)*mH*mCosDipAngle; +inline double Helix::cy(double s) const { + if (mSingularity) + return mCosDipAngle * mCosPhase; + else + return cos(mPhase + s * mH * mCurvature * mCosDipAngle) * mH * mCosDipAngle; } -inline double Helix::cz(double /* s */) const -{ - return mSinDipAngle; -} +inline double Helix::cz(double /* s */) const { return mSinDipAngle; } -inline const edm4hep::Vector3f& Helix::origin() const {return mOrigin;} +inline const edm4hep::Vector3f& Helix::origin() const { return mOrigin; } -inline edm4hep::Vector3f Helix::at(double s) const -{ - return edm4hep::Vector3f(x(s), y(s), z(s)); -} +inline edm4hep::Vector3f Helix::at(double s) const { return edm4hep::Vector3f(x(s), y(s), z(s)); } -inline edm4hep::Vector3f Helix::cat(double s) const -{ - return edm4hep::Vector3f(cx(s), cy(s), cz(s)); +inline edm4hep::Vector3f Helix::cat(double s) const { + return edm4hep::Vector3f(cx(s), cy(s), cz(s)); } -inline double Helix::pathLength(double X, double Y) const -{ - return fudgePathLength(edm4hep::Vector3f(X, Y, 0)); +inline double Helix::pathLength(double X, double Y) const { + return fudgePathLength(edm4hep::Vector3f(X, Y, 0)); } -inline int Helix::bad(double WorldSize) const -{ +inline int Helix::bad(double WorldSize) const { -// int ierr; - if (!::finite(mDipAngle )) return 11; - if (!::finite(mCurvature )) return 12; + // int ierr; + if (!::finite(mDipAngle)) + return 11; + if (!::finite(mCurvature)) + return 12; -// ierr = mOrigin.bad(WorldSize); -// if (ierr) return 3+ierr*100; + // ierr = mOrigin.bad(WorldSize); + // if (ierr) return 3+ierr*100; - if (::fabs(mDipAngle) >1.58) return 21; - double qwe = ::fabs(::fabs(mDipAngle)-M_PI/2); - if (qwe < 1./WorldSize ) return 31; + if (::fabs(mDipAngle) > 1.58) + return 21; + double qwe = ::fabs(::fabs(mDipAngle) - M_PI / 2); + if (qwe < 1. / WorldSize) + return 31; - if (::fabs(mCurvature) > WorldSize) return 22; - if (mCurvature < 0 ) return 32; + if (::fabs(mCurvature) > WorldSize) + return 22; + if (mCurvature < 0) + return 32; - if (abs(mH) != 1 ) return 24; + if (abs(mH) != 1) + return 24; - return 0; + return 0; } -} +} // namespace eicrecon From 57ccb508501b70310ad046a41a1ae4804da4e837 Mon Sep 17 00:00:00 2001 From: Xin Dong Date: Mon, 20 Oct 2025 19:17:53 -0400 Subject: [PATCH 12/17] clean up and remove unused variable --- src/algorithms/reco/Helix.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/algorithms/reco/Helix.cc b/src/algorithms/reco/Helix.cc index 5534bc31ca..242bb34229 100644 --- a/src/algorithms/reco/Helix.cc +++ b/src/algorithms/reco/Helix.cc @@ -372,7 +372,7 @@ double Helix::pathLength(const edm4hep::Vector3f& r, const edm4hep::Vector3f& n) double a, f, fp; double sOld = s = 0; - double shiftOld = 0; +// double shiftOld = 0; double shift; // (cos(angMax)-1)/angMax = 0.1 const double angMax = 0.21; @@ -400,7 +400,7 @@ double Helix::pathLength(const edm4hep::Vector3f& r, const edm4hep::Vector3f& n) shift = f / fp; } s -= shift; - shiftOld = shift; +// shiftOld = shift; if (fabs(sOld - s) < MaxPrecisionNeeded) break; sOld = s; From 69322bf6131899eb3d30b7bfc55f114ef98cc129 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 20 Oct 2025 23:18:21 +0000 Subject: [PATCH 13/17] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/algorithms/reco/Helix.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/algorithms/reco/Helix.cc b/src/algorithms/reco/Helix.cc index 242bb34229..575e71da29 100644 --- a/src/algorithms/reco/Helix.cc +++ b/src/algorithms/reco/Helix.cc @@ -372,7 +372,7 @@ double Helix::pathLength(const edm4hep::Vector3f& r, const edm4hep::Vector3f& n) double a, f, fp; double sOld = s = 0; -// double shiftOld = 0; + // double shiftOld = 0; double shift; // (cos(angMax)-1)/angMax = 0.1 const double angMax = 0.21; @@ -400,7 +400,7 @@ double Helix::pathLength(const edm4hep::Vector3f& r, const edm4hep::Vector3f& n) shift = f / fp; } s -= shift; -// shiftOld = shift; + // shiftOld = shift; if (fabs(sOld - s) < MaxPrecisionNeeded) break; sOld = s; From 7dcc3e0f49b1db51e3abc97e1ed1fcedf6f1b7c1 Mon Sep 17 00:00:00 2001 From: Xin Dong Date: Mon, 20 Oct 2025 19:32:42 -0400 Subject: [PATCH 14/17] further clean up --- src/algorithms/reco/Helix.h | 2 +- src/algorithms/reco/SecondaryVerticesHelix.cc | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/algorithms/reco/Helix.h b/src/algorithms/reco/Helix.h index 8e7f710c60..3631ad185c 100644 --- a/src/algorithms/reco/Helix.h +++ b/src/algorithms/reco/Helix.h @@ -107,7 +107,7 @@ class Helix { int bad(double world = 1.e+5) const; /// move the origin along the helix to s which becomes then s=0 - virtual void moveOrigin(double s); + void moveOrigin(double s); static const double NoSolution; diff --git a/src/algorithms/reco/SecondaryVerticesHelix.cc b/src/algorithms/reco/SecondaryVerticesHelix.cc index 40c65e4c8e..dd71037875 100644 --- a/src/algorithms/reco/SecondaryVerticesHelix.cc +++ b/src/algorithms/reco/SecondaryVerticesHelix.cc @@ -112,10 +112,10 @@ void SecondaryVerticesHelix::process(const SecondaryVerticesHelix::Input& input, double e2 = std::hypot(edm4hep::utils::magnitude(h2MomAtDca), particleSvc.particle(p2.getPDG()).mass); double pairE = e1 + e2; - double pairP = edm4hep::utils::magnitude(pairMom); +// double pairP = edm4hep::utils::magnitude(pairMom); - double m_inv2 = pairE * pairE - pairP * pairP; - double m_inv = (m_inv2 > 0) ? sqrt(m_inv2) : 0.; +// double m_inv2 = pairE * pairE - pairP * pairP; +// double m_inv = (m_inv2 > 0) ? sqrt(m_inv2) : 0.; double angle = edm4hep::utils::angleBetween(pairMom, pairPos - pVtxPos); if (cos(angle) < m_cfg.minCostheta) continue; From 71232c61d5c8a3cf2a8e27a06584e5cf9e7ef2f3 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 20 Oct 2025 23:33:13 +0000 Subject: [PATCH 15/17] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/algorithms/reco/SecondaryVerticesHelix.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/algorithms/reco/SecondaryVerticesHelix.cc b/src/algorithms/reco/SecondaryVerticesHelix.cc index dd71037875..f4fefe6b34 100644 --- a/src/algorithms/reco/SecondaryVerticesHelix.cc +++ b/src/algorithms/reco/SecondaryVerticesHelix.cc @@ -112,11 +112,11 @@ void SecondaryVerticesHelix::process(const SecondaryVerticesHelix::Input& input, double e2 = std::hypot(edm4hep::utils::magnitude(h2MomAtDca), particleSvc.particle(p2.getPDG()).mass); double pairE = e1 + e2; -// double pairP = edm4hep::utils::magnitude(pairMom); + // double pairP = edm4hep::utils::magnitude(pairMom); -// double m_inv2 = pairE * pairE - pairP * pairP; -// double m_inv = (m_inv2 > 0) ? sqrt(m_inv2) : 0.; - double angle = edm4hep::utils::angleBetween(pairMom, pairPos - pVtxPos); + // double m_inv2 = pairE * pairE - pairP * pairP; + // double m_inv = (m_inv2 > 0) ? sqrt(m_inv2) : 0.; + double angle = edm4hep::utils::angleBetween(pairMom, pairPos - pVtxPos); if (cos(angle) < m_cfg.minCostheta) continue; From 488a05c09adc4b395490d841656ae67f44dac3e7 Mon Sep 17 00:00:00 2001 From: Xin Dong Date: Mon, 20 Oct 2025 20:06:31 -0400 Subject: [PATCH 16/17] added header guard --- src/algorithms/reco/Helix.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/algorithms/reco/Helix.h b/src/algorithms/reco/Helix.h index 3631ad185c..4080394f1c 100644 --- a/src/algorithms/reco/Helix.h +++ b/src/algorithms/reco/Helix.h @@ -1,6 +1,8 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) 2024 Xin Dong, Rongrong Ma +#pragma once + #include #include #include From dbc662f2b9377608e8713da02965be54625a5279 Mon Sep 17 00:00:00 2001 From: Wouter Deconinck Date: Mon, 20 Oct 2025 19:46:17 -0500 Subject: [PATCH 17/17] SecondaryVertices factory based on Helix method (fix: iwyu) (#2145) This PR applies the include-what-you-use fixes as suggested by https://github.com/eic/EICrecon/actions/runs/18668462362. Please merge this PR into the branch `pr/secondaryvertex-helix` to resolve failures in PR #2144. Auto-generated by [create-pull-request][1] [1]: https://github.com/peter-evans/create-pull-request Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/algorithms/reco/Helix.cc | 22 +++++++++---------- src/algorithms/reco/Helix.h | 20 +++++------------ src/algorithms/reco/SecondaryVerticesHelix.cc | 10 ++++----- 3 files changed, 20 insertions(+), 32 deletions(-) diff --git a/src/algorithms/reco/Helix.cc b/src/algorithms/reco/Helix.cc index 575e71da29..c59e421230 100644 --- a/src/algorithms/reco/Helix.cc +++ b/src/algorithms/reco/Helix.cc @@ -1,24 +1,22 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) 2024 Xin Dong, Rongrong Ma +#include +#include +#include +#include +#include +#include +#include +#include #include #include -#include +#include #include -#include -#include #include -#include -#include -#include -#include +#include -#include -#include #include "algorithms/reco/Helix.h" -#include -#include -#include namespace eicrecon { diff --git a/src/algorithms/reco/Helix.h b/src/algorithms/reco/Helix.h index 4080394f1c..ac52e49ce5 100644 --- a/src/algorithms/reco/Helix.h +++ b/src/algorithms/reco/Helix.h @@ -3,22 +3,14 @@ #pragma once -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include #include -#include #include +#include +#include +#include +#include +#include +#include namespace eicrecon { diff --git a/src/algorithms/reco/SecondaryVerticesHelix.cc b/src/algorithms/reco/SecondaryVerticesHelix.cc index f4fefe6b34..2f162fc3bf 100644 --- a/src/algorithms/reco/SecondaryVerticesHelix.cc +++ b/src/algorithms/reco/SecondaryVerticesHelix.cc @@ -1,22 +1,20 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) 2024 Daniel Brandenburg, Xin Dong -#include -#include +#include #include #include +#include #include +#include #include -#include #include -#include #include -#include #include #include -#include "algorithms/reco/SecondaryVerticesHelix.h" #include "algorithms/reco/Helix.h" +#include "algorithms/reco/SecondaryVerticesHelix.h" #include "algorithms/reco/SecondaryVerticesHelixConfig.h" #include "services/particle/ParticleSvc.h"