SecondaryVertices factory based on Helix method

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 <edm4eic/ReconstructedParticleCollection.h>
+#include <edm4eic/TrackParametersCollection.h>
+#include <edm4eic/unit_system.h>
+#include <edm4hep/Vector3f.h>
+#include <stdlib.h>
+#include <cmath>
+#include <iterator>
+#include <utility>
+
+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<double, double> pathLength(double r) const;
+
+  /// path length at given r (cylindrical r, cylinder axis at x,y)
+  std::pair<double, double> 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<double, double> 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

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 <Evaluator/DD4hepUnits.h>
+#include <edm4eic/VertexCollection.h>
+#include <edm4eic/unit_system.h>
+#include <edm4hep/Vector3f.h>
+#include <edm4hep/Vector4f.h>
+#include <edm4hep/utils/vector_utils.h>
+#include <fmt/core.h>
+#include <cmath>
+#include <gsl/pointers>
+#include <utility>
+#include <vector>
+
+#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<Helix> hVec;
+  hVec.clear();
+  std::vector<unsigned int> 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<double, double> 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