diff --git a/src/algorithms/reco/Helix.cc b/src/algorithms/reco/Helix.cc new file mode 100644 index 0000000000..c59e421230 --- /dev/null +++ b/src/algorithms/reco/Helix.cc @@ -0,0 +1,616 @@ +// 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 "algorithms/reco/Helix.h" + +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 < 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; + } + + // + // 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; + } + } + return s; +} + +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; + + // + // 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 < 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) { + // + // 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 < 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 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); + } +} +/* +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 diff --git a/src/algorithms/reco/Helix.h b/src/algorithms/reco/Helix.h new file mode 100644 index 0000000000..ac52e49ce5 --- /dev/null +++ b/src/algorithms/reco/Helix.h @@ -0,0 +1,222 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Xin Dong, Rongrong Ma + +#pragma once + +#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 + 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 + +// +// 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; +} + +} // namespace eicrecon diff --git a/src/algorithms/reco/SecondaryVerticesHelix.cc b/src/algorithms/reco/SecondaryVerticesHelix.cc new file mode 100644 index 0000000000..2f162fc3bf --- /dev/null +++ b/src/algorithms/reco/SecondaryVerticesHelix.cc @@ -0,0 +1,148 @@ +// 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 "algorithms/reco/Helix.h" +#include "algorithms/reco/SecondaryVerticesHelix.h" +#include "algorithms/reco/SecondaryVerticesHelixConfig.h" +#include "services/particle/ParticleSvc.h" + +namespace eicrecon { + +/** + * @brief Initialize the SecondaryVerticesHelix Algorithm + * + */ +void SecondaryVerticesHelix::init() {} + +/** + * @brief Produce a list of secondary vertex candidates + * + * @param rcvtx - input collection of all vertex candidates + * @return edm4eic::VertexCollection + */ +void SecondaryVerticesHelix::process(const SecondaryVerticesHelix::Input& input, + const SecondaryVerticesHelix::Output& output) const { + 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!"); + 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 hVec; + hVec.clear(); + std::vector indexVec; + indexVec.clear(); + for (unsigned int i = 0; const auto& p : *rcparts) { + 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; + } + + 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); + + } // end i2 + } // end i1 + +} // end process + +} // namespace eicrecon diff --git a/src/algorithms/reco/SecondaryVerticesHelix.h b/src/algorithms/reco/SecondaryVerticesHelix.h new file mode 100644 index 0000000000..cf065d2ec6 --- /dev/null +++ b/src/algorithms/reco/SecondaryVerticesHelix.h @@ -0,0 +1,38 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2025 Xin Dong + +#pragma once + +#include +#include +#include +#include // for basic_string +#include // for string_view + +#include "algorithms/interfaces/WithPodConfig.h" +#include "algorithms/reco/SecondaryVerticesHelixConfig.h" + +namespace eicrecon { + +using SecondaryVerticesHelixAlgorithm = algorithms::Algorithm< + algorithms::Input, + algorithms::Output>; + +class SecondaryVerticesHelix : public SecondaryVerticesHelixAlgorithm, + public WithPodConfig { + +public: + SecondaryVerticesHelix(std::string_view name) + : SecondaryVerticesHelixAlgorithm{ + name, + {"inputVertices", "inputParticles"}, + {"outputSecondaryVertices"}, + "Reconstruct secondary vertices in SecondaryVertices collection"} {} + + void init() final; + void process(const Input&, const Output&) const final; + +private: + SecondaryVerticesHelixConfig m_cfg; +}; +} // namespace eicrecon diff --git a/src/algorithms/reco/SecondaryVerticesHelixConfig.h b/src/algorithms/reco/SecondaryVerticesHelixConfig.h new file mode 100644 index 0000000000..cef2d8e300 --- /dev/null +++ b/src/algorithms/reco/SecondaryVerticesHelixConfig.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 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 minCostheta = 0.8; // costheta, theta: angle of V0 decay direction and momentum +}; + +} // namespace eicrecon diff --git a/src/factories/reco/SecondaryVerticesHelix_factory.h b/src/factories/reco/SecondaryVerticesHelix_factory.h new file mode 100644 index 0000000000..b52a459272 --- /dev/null +++ b/src/factories/reco/SecondaryVerticesHelix_factory.h @@ -0,0 +1,51 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Xin Dong + +#pragma once + +#include +#include +#include +#include +#include + +#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 SecondaryVerticesHelix_factory + : public JOmniFactory { + +public: + using AlgoT = eicrecon::SecondaryVerticesHelix; + +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 531ba981c3..d8469305d1 100644 --- a/src/global/reco/reco.cc +++ b/src/global/reco/reco.cc @@ -37,6 +37,7 @@ #include "factories/reco/MC2ReconstructedParticle_factory.h" #include "factories/reco/MatchClusters_factory.h" #include "factories/reco/PrimaryVertices_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" @@ -73,7 +74,6 @@ void InitPlugin(JApplication* app) { "EcalEndcapPClusterAssociations"}, {"EcalClusterAssociations"}, app)); - // Create ReconstructedParticles app->Add(new JOmniFactoryGeneratorT( "ReconstructedParticlesWithAssoc", { @@ -267,5 +267,9 @@ void InitPlugin(JApplication* app) { app->Add(new JOmniFactoryGeneratorT( "PrimaryVertices", {"CentralTrackVertices"}, {"PrimaryVertices"}, {}, 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 3bf936a201..1fa5c4d45e 100644 --- a/src/services/io/podio/JEventProcessorPODIO.cc +++ b/src/services/io/podio/JEventProcessorPODIO.cc @@ -261,6 +261,7 @@ JEventProcessorPODIO::JEventProcessorPODIO() { "ScatteredElectronsTruth", "ScatteredElectronsEMinusPz", "PrimaryVertices", + "SecondaryVerticesHelix", "BarrelClusters", "HadronicFinalState",