algorithms/reco/FarForwardLambdaReconstruction.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2025 Sebouh Paul
0003
0004 #include <Evaluator/DD4hepUnits.h>
0005 #include <TVector3.h>
0006 #include <edm4eic/ReconstructedParticleCollection.h>
0007 #include <edm4hep/Vector3f.h>
0008 #include <edm4hep/utils/vector_utils.h>
0009 #include <fmt/core.h>
0010 #include <cmath>
0011 #include <cstddef>
0012 #include <gsl/pointers>
0013 #include <stdexcept>
0014 #include <vector>
0015
0016 #include "FarForwardLambdaReconstruction.h"
0017 #include "TLorentzVector.h"
0018 /**
0019 Creates Lambda candidates from a neutron and two photons from a pi0 decay
0020  */
0021
0022 namespace eicrecon {
0023
0024 void FarForwardLambdaReconstruction::init() {
0025
0026   try {
0027     m_zMax = m_detector->constant<double>(m_cfg.offsetPositionName);
0028   } catch (std::runtime_error&) {
0029     m_zMax = 35800; // default value
0030     trace("Failed to get {} from the detector, using default value of {}", m_cfg.offsetPositionName,
0031           m_zMax);
0032   }
0033 }
0034
0035 /* converts one type of vector format to another */
0036 void toTVector3(TVector3& v1, const edm4hep::Vector3f& v2) { v1.SetXYZ(v2.x, v2.y, v2.z); }
0037
0038 void FarForwardLambdaReconstruction::process(
0039     const FarForwardLambdaReconstruction::Input& input,
0040     const FarForwardLambdaReconstruction::Output& output) const {
0041   const auto [neutrals]                                = input;
0042   auto [out_lambdas, out_decay_products]               = output;
0043   std::vector<edm4eic::ReconstructedParticle> neutrons = {};
0044   std::vector<edm4eic::ReconstructedParticle> gammas   = {};
0045   for (auto part : *neutrals) {
0046     if (part.getPDG() == 2112) {
0047       neutrons.push_back(part);
0048     }
0049     if (part.getPDG() == 22) {
0050       gammas.push_back(part);
0051     }
0052   }
0053
0054   if (neutrons.empty() || gammas.size() < 2) {
0055     return;
0056   }
0057
0058   static const double m_neutron = m_particleSvc.particle(2112).mass;
0059   static const double m_pi0     = m_particleSvc.particle(111).mass;
0060   static const double m_lambda  = m_particleSvc.particle(3122).mass;
0061
0062   for (std::size_t i_n = 0; i_n < neutrons.size(); i_n++) {
0063     for (std::size_t i_1 = 0; i_1 < gammas.size() - 1; i_1++) {
0064       for (std::size_t i_2 = i_1 + 1; i_2 < gammas.size(); i_2++) {
0065
0066         double En = neutrons[i_n].getEnergy();
0067         double pn = sqrt(En * En - m_neutron * m_neutron);
0068         double E1 = gammas[i_1].getEnergy();
0069         double E2 = gammas[i_2].getEnergy();
0070         TVector3 xn;
0071         TVector3 x1;
0072         TVector3 x2;
0073
0074         toTVector3(xn, neutrons[0].getReferencePoint() * dd4hep::mm);
0075         toTVector3(x1, gammas[0].getReferencePoint() * dd4hep::mm);
0076         toTVector3(x2, gammas[1].getReferencePoint() * dd4hep::mm);
0077
0078         xn.RotateY(-m_cfg.globalToProtonRotation);
0079         x1.RotateY(-m_cfg.globalToProtonRotation);
0080         x2.RotateY(-m_cfg.globalToProtonRotation);
0081
0082         debug("nx recon = {}, g1x recon = {}, g2x recon = {}", xn.X(), x1.X(), x2.X());
0083         debug("nz recon = {}, g1z recon = {}, g2z recon = {}, z face = {}", xn.Z(), x1.Z(), x2.Z(),
0084               m_zMax);
0085
0086         TVector3 vtx(0, 0, 0);
0087         double f                   = 0;
0088         double df                  = 0.5;
0089         double theta_open_expected = 2 * asin(m_pi0 / (2 * sqrt(E1 * E2)));
0090         TLorentzVector n;
0091         TLorentzVector g1;
0092         TLorentzVector g2;
0093         TLorentzVector lambda;
0094         for (int i = 0; i < m_cfg.iterations; i++) {
0095           n                 = {pn * (xn - vtx).Unit(), En};
0096           g1                = {E1 * (x1 - vtx).Unit(), E1};
0097           g2                = {E2 * (x2 - vtx).Unit(), E2};
0098           double theta_open = g1.Angle(g2.Vect());
0099           lambda            = n + g1 + g2;
0100           if (theta_open > theta_open_expected) {
0101             f -= df;
0102           } else if (theta_open < theta_open_expected) {
0103             f += df;
0104           }
0105
0106           vtx = lambda.Vect() * (f * m_zMax / lambda.Z());
0107           df /= 2;
0108         }
0109
0110         double mass_rec = lambda.M();
0111         if (std::abs(mass_rec - m_lambda) > m_cfg.lambdaMaxMassDev) {
0112           continue;
0113         }
0114
0115         // rotate everything back to the lab coordinates.
0116         vtx.RotateY(m_cfg.globalToProtonRotation);
0117         lambda.RotateY(m_cfg.globalToProtonRotation);
0118         n.RotateY(m_cfg.globalToProtonRotation);
0119         g1.RotateY(m_cfg.globalToProtonRotation);
0120         g2.RotateY(m_cfg.globalToProtonRotation);
0121
0122         auto b = -lambda.BoostVector();
0123         n.Boost(b);
0124         g1.Boost(b);
0125         g2.Boost(b);
0126
0127         //convert vertex to mm:
0128         vtx = vtx * (1 / dd4hep::mm);
0129
0130         auto rec_lambda = out_lambdas->create();
0131         rec_lambda.setPDG(3122);
0132
0133         rec_lambda.setEnergy(lambda.E());
0134         rec_lambda.setMomentum({static_cast<float>(lambda.X()), static_cast<float>(lambda.Y()),
0135                                 static_cast<float>(lambda.Z())});
0136         rec_lambda.setReferencePoint({static_cast<float>(vtx.X()), static_cast<float>(vtx.Y()),
0137                                       static_cast<float>(vtx.Z())});
0138         rec_lambda.setCharge(0);
0139         rec_lambda.setMass(mass_rec);
0140
0141         auto neutron_cm = out_decay_products->create();
0142         neutron_cm.setPDG(2112);
0143         neutron_cm.setEnergy(n.E());
0144         neutron_cm.setMomentum(
0145             {static_cast<float>(n.X()), static_cast<float>(n.Y()), static_cast<float>(n.Z())});
0146         neutron_cm.setReferencePoint({static_cast<float>(vtx.X()), static_cast<float>(vtx.Y()),
0147                                       static_cast<float>(vtx.Z())});
0148         neutron_cm.setCharge(0);
0149         neutron_cm.setMass(m_neutron);
0150         //link the reconstructed lambda to the input neutron,
0151         // and the cm neutron to the input neutron
0152         rec_lambda.addToParticles(neutrons[i_n]);
0153         neutron_cm.addToParticles(neutrons[i_n]);
0154
0155         auto gamma1_cm = out_decay_products->create();
0156         gamma1_cm.setPDG(22);
0157         gamma1_cm.setEnergy(g1.E());
0158         gamma1_cm.setMomentum(
0159             {static_cast<float>(g1.X()), static_cast<float>(g1.Y()), static_cast<float>(g1.Z())});
0160         gamma1_cm.setReferencePoint({static_cast<float>(vtx.X()), static_cast<float>(vtx.Y()),
0161                                      static_cast<float>(vtx.Z())});
0162         gamma1_cm.setCharge(0);
0163         gamma1_cm.setMass(0);
0164         rec_lambda.addToParticles(gammas[i_1]);
0165         gamma1_cm.addToParticles(gammas[i_1]);
0166
0167         auto gamma2_cm = out_decay_products->create();
0168         gamma2_cm.setPDG(22);
0169         gamma2_cm.setEnergy(g2.E());
0170         gamma2_cm.setMomentum(
0171             {static_cast<float>(g2.X()), static_cast<float>(g2.Y()), static_cast<float>(g2.Z())});
0172         gamma2_cm.setReferencePoint({static_cast<float>(vtx.X()), static_cast<float>(vtx.Y()),
0173                                      static_cast<float>(vtx.Z())});
0174         gamma2_cm.setCharge(0);
0175         gamma2_cm.setMass(0);
0176         rec_lambda.addToParticles(gammas[i_2]);
0177         gamma2_cm.addToParticles(gammas[i_2]);
0178       }
0179     }
0180   }
0181 }
0182 } // namespace eicrecon