algorithms/reco/SecondaryVerticesHelix.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2025 Daniel Brandenburg, Xin Dong
0003
0004 #include <DD4hep/Fields.h>
0005 #include <Evaluator/DD4hepUnits.h>
0006 #include <Math/GenVector/Cartesian3D.h>
0007 #include <Math/GenVector/DisplacementVector3D.h>
0008 #include <edm4eic/VertexCollection.h>
0009 #include <edm4eic/unit_system.h>
0010 #include <edm4hep/Vector3f.h>
0011 #include <edm4hep/Vector4f.h>
0012 #include <edm4hep/utils/vector_utils.h>
0013 #include <fmt/core.h>
0014 #include <cmath>
0015 #include <gsl/pointers>
0016 #include <utility>
0017 #include <vector>
0018
0019 #include "algorithms/reco/Helix.h"
0020 #include "algorithms/reco/SecondaryVerticesHelix.h"
0021 #include "algorithms/reco/SecondaryVerticesHelixConfig.h"
0022 #include "services/particle/ParticleSvc.h"
0023
0024 namespace eicrecon {
0025
0026 /**
0027    * @brief Initialize the SecondaryVerticesHelix Algorithm
0028    *
0029    */
0030 void SecondaryVerticesHelix::init() {}
0031
0032 /**
0033    * @brief Produce a list of secondary vertex candidates
0034    *
0035    * @param rcvtx  - input collection of all vertex candidates
0036    * @return edm4eic::VertexCollection
0037    */
0038 void SecondaryVerticesHelix::process(const SecondaryVerticesHelix::Input& input,
0039                                      const SecondaryVerticesHelix::Output& output) const {
0040   const auto [rcvtx, rcparts]   = input;
0041   auto [out_secondary_vertices] = output;
0042
0043   auto& particleSvc = algorithms::ParticleSvc::instance();
0044
0045   if ((*rcvtx).size() == 0) {
0046     info("No primary vertex in this event! Skip secondary vertex finder!");
0047     return;
0048   }
0049   const auto pVtxPos4f = (*rcvtx)[0].getPosition();
0050   // convert to cm
0051   edm4hep::Vector3f pVtxPos(pVtxPos4f.x * edm4eic::unit::mm / edm4eic::unit::cm,
0052                             pVtxPos4f.y * edm4eic::unit::mm / edm4eic::unit::cm,
0053                             pVtxPos4f.z * edm4eic::unit::mm / edm4eic::unit::cm);
0054
0055   auto fieldObj = m_det->field();
0056   auto field    = fieldObj.magneticField(
0057       {pVtxPos4f.x / edm4eic::unit::mm * dd4hep::mm, pVtxPos4f.y / edm4eic::unit::mm * dd4hep::mm,
0058           pVtxPos4f.z / edm4eic::unit::mm * dd4hep::mm}); // in unit of dd4hep::tesla
0059   float b_field = field.z();
0060
0061   info("\tPrimary vertex = ({},{},{})cm \t b field = {} tesla", pVtxPos.x, pVtxPos.y, pVtxPos.z,
0062        b_field / dd4hep::tesla);
0063
0064   std::vector<Helix> hVec;
0065   hVec.clear();
0066   std::vector<unsigned int> indexVec;
0067   indexVec.clear();
0068   for (unsigned int i = 0; const auto& p : *rcparts) {
0069     if (p.getCharge() == 0)
0070       continue;
0071     Helix h(p, b_field);
0072     double dca = h.distance(pVtxPos) * edm4eic::unit::cm;
0073     if (dca < m_cfg.minDca)
0074       continue;
0075
0076     hVec.push_back(h);
0077     indexVec.push_back(i);
0078     ++i;
0079   }
0080
0081   if (hVec.size() != indexVec.size())
0082     return;
0083
0084   debug("\tVector size {}, {}", hVec.size(), indexVec.size());
0085
0086   for (unsigned int i1 = 0; i1 < hVec.size(); ++i1) {
0087     for (unsigned int i2 = i1 + 1; i2 < hVec.size(); ++i2) {
0088       const auto& p1 = (*rcparts)[indexVec[i1]];
0089       const auto& p2 = (*rcparts)[indexVec[i2]];
0090
0091       if (!(m_cfg.unlikesign && p1.getCharge() + p2.getCharge() == 0))
0092         continue;
0093
0094       const auto& h1 = hVec[i1];
0095       const auto& h2 = hVec[i2];
0096
0097       // Helix function uses cm unit
0098       double dca1 = h1.distance(pVtxPos) * edm4eic::unit::cm;
0099       double dca2 = h2.distance(pVtxPos) * edm4eic::unit::cm;
0100       if (dca1 < m_cfg.minDca || dca2 < m_cfg.minDca)
0101         continue;
0102
0103       std::pair<double, double> const ss = h1.pathLengths(h2);
0104       edm4hep::Vector3f h1AtDcaTo2       = h1.at(ss.first);
0105       edm4hep::Vector3f h2AtDcaTo1       = h2.at(ss.second);
0106
0107       double dca12 = edm4hep::utils::magnitude(h1AtDcaTo2 - h2AtDcaTo1) * edm4eic::unit::cm;
0108       if (std::isnan(dca12))
0109         continue;
0110       if (dca12 > m_cfg.maxDca12)
0111         continue;
0112       edm4hep::Vector3f pairPos = 0.5 * (h1AtDcaTo2 + h2AtDcaTo1);
0113
0114       edm4hep::Vector3f h1MomAtDca = h1.momentumAt(ss.first, b_field);
0115       edm4hep::Vector3f h2MomAtDca = h2.momentumAt(ss.second, b_field);
0116       edm4hep::Vector3f pairMom    = h1MomAtDca + h2MomAtDca;
0117
0118       double e1 =
0119           std::hypot(edm4hep::utils::magnitude(h1MomAtDca), particleSvc.particle(p1.getPDG()).mass);
0120       double e2 =
0121           std::hypot(edm4hep::utils::magnitude(h2MomAtDca), particleSvc.particle(p2.getPDG()).mass);
0122       double pairE = e1 + e2;
0123       double angle = edm4hep::utils::angleBetween(pairMom, pairPos - pVtxPos);
0124       if (cos(angle) < m_cfg.minCostheta)
0125         continue;
0126
0127       double beta = edm4hep::utils::magnitude(pairMom) / pairE;
0128       double time = edm4hep::utils::magnitude(pairPos - pVtxPos) / (beta * dd4hep::c_light);
0129       edm4hep::Vector3f dL = pairPos - pVtxPos; // in cm
0130       edm4hep::Vector3f decayL(dL.x * edm4eic::unit::cm, dL.y * edm4eic::unit::cm,
0131                                dL.z * edm4eic::unit::cm);
0132       double dca2pv = edm4hep::utils::magnitude(decayL) * sin(angle);
0133       if (dca2pv > m_cfg.maxDca)
0134         continue;
0135
0136       auto v0 = out_secondary_vertices->create();
0137       v0.setType(2); // 2 for secondary
0138       v0.setPosition({(float)(pairPos.x * edm4eic::unit::cm / edm4eic::unit::mm),
0139                       (float)(pairPos.y * edm4eic::unit::cm / edm4eic::unit::mm),
0140                       (float)(pairPos.z * edm4eic::unit::cm / edm4eic::unit::mm), (float)time});
0141       v0.addToAssociatedParticles(p1);
0142       v0.addToAssociatedParticles(p2);
0143
0144       info("One secondary vertex found at (x,y,z) = ({}, {}, {}) mm.",
0145            pairPos.x * edm4eic::unit::cm / edm4eic::unit::mm,
0146            pairPos.y * edm4eic::unit::cm / edm4eic::unit::mm,
0147            pairPos.x * edm4eic::unit::cm / edm4eic::unit::mm);
0148
0149     } // end i2
0150   } // end i1
0151
0152 } // end process
0153
0154 } // namespace eicrecon