algorithms/reco/FarForwardNeutralsReconstruction.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2026 Sebouh Paul, Baptiste Fraisse
0003
0004 #include <Evaluator/DD4hepUnits.h>
0005 #include <edm4eic/ClusterCollection.h>
0006 #include <edm4eic/ReconstructedParticleCollection.h>
0007 #include <edm4hep/Vector3f.h>
0008 #include <edm4hep/utils/vector_utils.h>
0009 #include <stddef.h>
0010 #include <algorithm>
0011 #include <cmath>
0012 #include <stdexcept>
0013 #include <tuple>
0014 #include <vector>
0015
0016 #include "FarForwardNeutralsReconstruction.h"
0017
0018 namespace eicrecon {
0019
0020 void FarForwardNeutralsReconstruction::init() {
0021   try {
0022     m_gammaZMax =
0023         m_cfg.gammaZMaxOffset + m_detector->constant<double>(m_cfg.offsetPositionName) / dd4hep::mm;
0024   } catch (std::runtime_error&) {
0025     m_gammaZMax = m_cfg.gammaZMaxOffset + 35800.0;
0026     trace("Failed to get {} from the detector, using default value of {}", m_cfg.offsetPositionName,
0027           m_gammaZMax);
0028   }
0029   trace("gamma detection params:   max length={},   max width={},   max z={}", m_cfg.gammaMaxLength,
0030         m_cfg.gammaMaxWidth, m_gammaZMax);
0031 }
0032
0033 /*
0034      check that the cluster position is within the correct range into ZDC,
0035      and that the sqrt(largest eigenvalue) is less than gamma_max_length,
0036      and that the sqrt(second largest eigenvalue) is less than gamma_max_width
0037 */
0038 bool FarForwardNeutralsReconstruction::isGamma(const edm4eic::Cluster& cluster) const {
0039
0040   double l1 = sqrt(cluster.getShapeParameters(4)) * dd4hep::mm;
0041   double l2 = sqrt(cluster.getShapeParameters(5)) * dd4hep::mm;
0042   double l3 = sqrt(cluster.getShapeParameters(6)) * dd4hep::mm;
0043
0044   // z in the local coordinates
0045   double z = (cluster.getPosition().z * cos(m_cfg.globalToProtonRotation) +
0046               cluster.getPosition().x * sin(m_cfg.globalToProtonRotation)) *
0047              dd4hep::mm;
0048
0049   trace("z recon = {}", z);
0050   trace("l1 = {}, l2 = {}, l3 = {}", l1, l2, l3);
0051
0052   bool isZMoreThanMax = (z > m_gammaZMax);
0053   bool isLengthMoreThanMax =
0054       (l1 > m_cfg.gammaMaxLength || l2 > m_cfg.gammaMaxLength || l3 > m_cfg.gammaMaxLength);
0055   bool areWidthsMoreThanMax = static_cast<int>(l1 > m_cfg.gammaMaxWidth) +
0056                                   static_cast<int>(l2 > m_cfg.gammaMaxWidth) +
0057                                   static_cast<int>(l3 > m_cfg.gammaMaxWidth) >=
0058                               2;
0059
0060   return !(isZMoreThanMax || isLengthMoreThanMax || areWidthsMoreThanMax);
0061 }
0062
0063 double FarForwardNeutralsReconstruction::corrPower(double E, const std::vector<double>& coeffs) {
0064
0065   if (coeffs.size() != 2) {
0066     throw std::runtime_error(
0067         "Energy correction for Lambda reconstruction requires 2 coefficients: a * E^b.");
0068   }
0069
0070   return coeffs[0] * std::pow(E, coeffs[1]);
0071 }
0072
0073 int FarForwardNeutralsReconstruction::processNeutralCalo(
0074     const edm4eic::ClusterCollection* clusters,
0075     edm4eic::ReconstructedParticleCollection* out_neutrals,
0076     const std::vector<double>& gammaScaleCoeff, const std::vector<double>& neutronScaleCoeff,
0077     bool canDetectGammas, bool canDetectNeutrons, const CorrFunc& gammaCorr,
0078     const CorrFunc& neutronCorr, GammaMode gammaMode, double gammaLeaderFracMin, double clusterEmin,
0079     NeutronMode neutronMode, bool associateAllClustersToNeutron) const {
0080
0081   (void)gammaLeaderFracMin;
0082
0083   const double m_neutron = m_particleSvc.particle(2112).mass;
0084
0085   if ((clusters == nullptr) || clusters->empty()) {
0086     return 0;
0087   }
0088   if (!canDetectGammas) {
0089     gammaMode = GammaMode::None;
0090   }
0091
0092   // gammas from clusters
0093   auto makeGamma = [&](const edm4eic::Cluster& cl) {
0094     auto rec = out_neutrals->create();
0095     rec.setPDG(22);
0096
0097     const auto pos = cl.getPosition();
0098     const double E = gammaCorr(cl.getEnergy(), gammaScaleCoeff);
0099
0100     const double r = edm4hep::utils::magnitude(pos);
0101     if (r > 0) {
0102       rec.setMomentum(pos * (E / r));
0103     }
0104
0105     rec.setEnergy(E);
0106     rec.setCharge(0);
0107     rec.setMass(0);
0108     rec.addToClusters(cl);
0109   };
0110
0111   std::vector<const edm4eic::Cluster*> gamma_used;
0112   gamma_used.reserve(clusters->size());
0113
0114   // gammaMode == LeaderOnly
0115   if (gammaMode == GammaMode::LeaderOnly) {
0116
0117     std::vector<int> idx;
0118     idx.reserve(clusters->size());
0119
0120     double Esum = 0.0;
0121     for (int i = 0, n = (int)clusters->size(); i < n; ++i) {
0122       const double E = (*clusters)[i].getEnergy();
0123       if (E < clusterEmin) {
0124         continue;
0125       }
0126       Esum += E;
0127       idx.push_back(i);
0128     }
0129
0130     if (idx.empty() || Esum <= 0.0) {
0131       return 0;
0132     }
0133
0134     // Sort clusters by decreasing energy before keeping the leading ones.
0135     std::ranges::sort(
0136         idx, [&](int a, int b) { return (*clusters)[a].getEnergy() > (*clusters)[b].getEnergy(); });
0137
0138     // Keep only the leading clusters in order to limit combinatorial
0139     const size_t Nkeep = 4;
0140     for (size_t k = 0; k < std::min(Nkeep, idx.size()); ++k) {
0141       const auto& cl = (*clusters)[idx[k]];
0142       makeGamma(cl);
0143       gamma_used.push_back(&cl);
0144     }
0145
0146   }
0147
0148   // gammaMode == AllPassing
0149   else if (gammaMode == GammaMode::AllPassing) {
0150     for (const auto& cl : *clusters) {
0151       const double E = cl.getEnergy();
0152       if (E < clusterEmin)
0153         continue;
0154       if (isGamma(cl)) {
0155         makeGamma(cl);
0156         gamma_used.push_back(&cl);
0157       }
0158     }
0159   }
0160
0161   // gammaMode == None => nothing
0162   auto is_used_as_gamma = [&](const edm4eic::Cluster& cl) {
0163     for (auto* p : gamma_used)
0164       if (p == &cl)
0165         return true;
0166     return false;
0167   };
0168
0169   // neutrons from clusters
0170   const edm4eic::Cluster* leaderN = nullptr;
0171   double E_leader                 = -1.0;
0172
0173   double E_sum = 0.0;
0174   std::vector<const edm4eic::Cluster*> kept;
0175   kept.reserve(clusters->size());
0176
0177   for (const auto& cl : *clusters) {
0178     if (gammaMode != GammaMode::None && is_used_as_gamma(cl))
0179       continue;
0180
0181     const double E = cl.getEnergy();
0182     if (E < clusterEmin)
0183       continue;
0184
0185     E_sum += E;
0186     kept.push_back(&cl);
0187
0188     if (E > E_leader) {
0189       E_leader = E;
0190       leaderN  = &cl;
0191     }
0192   }
0193
0194   if (!canDetectNeutrons)
0195     return 0;
0196
0197   double En_raw = 0.0;
0198   edm4hep::Vector3f n_pos{0, 0, 0};
0199
0200   if (neutronMode == NeutronMode::LeaderOnly) {
0201     if (!leaderN || E_leader <= 0.0)
0202       return 0;
0203     En_raw = E_leader;
0204     n_pos  = leaderN->getPosition();
0205
0206     kept.clear();
0207     kept.push_back(leaderN);
0208   } else if (neutronMode == NeutronMode::SumAll) {
0209     if (E_sum <= 0.0 || kept.empty())
0210       return 0;
0211     En_raw = E_sum;
0212     if (!leaderN || E_leader <= 0.0)
0213       return 0;
0214     n_pos = leaderN->getPosition();
0215   } else {
0216     return 0;
0217   }
0218
0219   // apply calibration laws
0220   const double En = neutronCorr(En_raw, neutronScaleCoeff);
0221
0222   auto rec = out_neutrals->create();
0223   rec.setPDG(2112);
0224   rec.setEnergy(En);
0225   rec.setCharge(0);
0226   rec.setMass(m_neutron);
0227
0228   const double r = edm4hep::utils::magnitude(n_pos);
0229   if (r > 0) {
0230     const double psq = std::max(0.0, En * En - m_neutron * m_neutron);
0231     rec.setMomentum(n_pos * (std::sqrt(psq) / r));
0232   }
0233   rec.setReferencePoint(n_pos);
0234
0235   if (associateAllClustersToNeutron) {
0236     for (const auto& cl : *clusters)
0237       rec.addToClusters(cl);
0238   } else {
0239     for (auto* clp : kept)
0240       rec.addToClusters(*clp);
0241   }
0242
0243   return 1;
0244 }
0245
0246 void FarForwardNeutralsReconstruction::process(
0247     const FarForwardNeutralsReconstruction::Input& input,
0248     const FarForwardNeutralsReconstruction::Output& output) const {
0249
0250   // Unpacking
0251   const auto [clustersHcal, clustersB0, clustersEcalEndcapP, clustersLFHCAL]           = input;
0252   auto [out_neutralsHcal, out_neutralsB0, out_neutralsEcalEndcapP, out_neutralsLFHCAL] = output;
0253
0254   // Global
0255   int n_neutrons = 0;
0256
0257   // ZDC-Hcal
0258   n_neutrons += processNeutralCalo(clustersHcal, out_neutralsHcal,
0259                                    /*gammaScaleCoeff=*/m_cfg.gammaScaleCorrCoeffHcalZDC,
0260                                    /*neutronScaleCoeff=*/m_cfg.neutronScaleCorrCoeffHcalZDC,
0261                                    /*canDetectGammas=*/true,
0262                                    /*canDetectNeutrons=*/true,
0263                                    /*gammaCorr=*/corrPower,
0264                                    /*neutronCorr=*/corrPower,
0265                                    /*gammaMode=*/GammaMode::AllPassing,
0266                                    /*gammaLeaderFracMin=*/0.0,
0267                                    /*clusterEmin=*/m_cfg.clusterEminHcalZDC,
0268                                    /*neutronMode=*/NeutronMode::SumAll,
0269                                    /*associateAllClustersToNeutron=*/true);
0270
0271   // B0-Ecal
0272   n_neutrons += processNeutralCalo(clustersB0, out_neutralsB0,
0273                                    /*gammaScaleCoeff=*/m_cfg.gammaScaleCorrCoeffB0Ecal,
0274                                    /*neutronScaleCoeff=*/m_cfg.neutronScaleCorrCoeffB0Ecal,
0275                                    /*canDetectGammas=*/true,
0276                                    /*canDetectNeutrons=*/false,
0277                                    /*gammaCorr=*/corrPower,
0278                                    /*neutronCorr=*/corrPower,
0279                                    /*gammaMode=*/GammaMode::LeaderOnly,
0280                                    /*gammaLeaderFracMin=*/0.0,
0281                                    /*clusterEmin=*/m_cfg.clusterEminB0Ecal,
0282                                    /*neutronMode=*/NeutronMode::None,
0283                                    /*associateAllClustersToNeutron=*/false);
0284
0285   // EndcapP-Ecal
0286   n_neutrons += processNeutralCalo(clustersEcalEndcapP, out_neutralsEcalEndcapP,
0287                                    /*gammaScaleCoeff=*/m_cfg.gammaScaleCorrCoeffEcalEndcapP,
0288                                    /*neutronScaleCoeff=*/m_cfg.neutronScaleCorrCoeffEcalEndcapP,
0289                                    /*canDetectGammas=*/true,
0290                                    /*canDetectNeutrons=*/false,
0291                                    /*gammaCorr=*/corrPower,
0292                                    /*neutronCorr=*/corrPower,
0293                                    /*gammaMode=*/GammaMode::LeaderOnly,
0294                                    /*gammaLeaderFracMin=*/0.0,
0295                                    /*clusterEmin=*/m_cfg.clusterEminEcalEndcapP,
0296                                    /*neutronMode=*/NeutronMode::None,
0297                                    /*associateAllClustersToNeutron=*/false);
0298
0299   // LFHCAL
0300   n_neutrons += processNeutralCalo(clustersLFHCAL, out_neutralsLFHCAL,
0301                                    /*gammaScaleCoeff=*/m_cfg.gammaScaleCorrCoeffLFHCAL,
0302                                    /*neutronScaleCoeff=*/m_cfg.neutronScaleCorrCoeffLFHCAL,
0303                                    /*canDetectGammas=*/false,
0304                                    /*canDetectNeutrons=*/true,
0305                                    /*gammaCorr=*/corrPower,
0306                                    /*neutronCorr=*/corrPower,
0307                                    /*gammaMode=*/GammaMode::None,
0308                                    /*gammaLeaderFracMin=*/0.0,
0309                                    /*clusterEmin=*/m_cfg.clusterEminLFHCAL,
0310                                    /*neutronMode=*/NeutronMode::LeaderOnly,
0311                                    /*associateAllClustersToNeutron=*/false);
0312
0313   debug("Found {} neutron candidates", n_neutrons);
0314 }
0315
0316 } // namespace eicrecon