algorithms/reco/FarForwardNeutralsReconstruction.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2025 Sebouh Paul
0003 #include <Evaluator/DD4hepUnits.h>
0004 #include <edm4eic/ClusterCollection.h>
0005 #include <edm4eic/ReconstructedParticleCollection.h>
0006 #include <edm4hep/Vector3f.h>
0007 #include <edm4hep/utils/vector_utils.h>
0008 #include <fmt/core.h>
0009 #include <cmath>
0010 #include <gsl/pointers>
0011 #include <stdexcept>
0012 #include <vector>
0013
0014 #include "FarForwardNeutralsReconstruction.h"
0015
0016 /**
0017  Creates photon candidate Reconstructed Particles, using clusters which fulfill cuts on the position of their CoG position, length (sqrt of largest eigenvector of their moment matrix), and width (sqrt of second largest eigenvector of their moment matrix).  Its energy is the energy of the cluster times a correction factor.
0018  Also creates a "neutron candidate" Reconstructed Particle consisting of all remaining clusters in the collection. Its energy is the sum of the energies of the constituent clusters
0019  times a correction factor, and its direction is the direction from the origin to the position
0020  of the most energetic cluster.  The correction factors in both cases are given by 1/(1+c[0]+c[1]/sqrt(E)+c[2]/E),
0021  where c is the coefficients and E is the uncorrected energy in GeV.  Different correction coefficients are used for photons vs for neutrons.  This form was chosen
0022  empirically based on the discrepancies in single-neutron and single-photon MC simulations between the uncorrected
0023  reconstructed energies and the truth energies of the neutrons.  The parameter values should be co-tuned with those of the clustering algorithm being used.
0024  */
0025
0026 namespace eicrecon {
0027
0028 void FarForwardNeutralsReconstruction::init() {
0029
0030   try {
0031     m_gammaZMax =
0032         m_cfg.gammaZMaxOffset + m_detector->constant<double>(m_cfg.offsetPositionName) / dd4hep::mm;
0033   } catch (std::runtime_error&) {
0034     m_gammaZMax = m_cfg.gammaZMaxOffset + 35800;
0035     trace("Failed to get {} from the detector, using default value of {}", m_cfg.offsetPositionName,
0036           m_gammaZMax);
0037   }
0038
0039   if (m_cfg.neutronScaleCorrCoeffHcal.size() < 3) {
0040     error("Invalid configuration.  m_cfg.neutronScaleCorrCoeffHcal should have at least 3 "
0041           "parameters");
0042     throw std::runtime_error("Invalid configuration.  m_cfg.neutronScaleCorrCoeffHcal should have "
0043                              "at least 3 parameters");
0044   }
0045   if (m_cfg.gammaScaleCorrCoeffHcal.size() < 3) {
0046     error("Invalid configuration.  m_cfg.gamma_scale_corr_coeff_ecal should have at least 3 "
0047           "parameters");
0048     throw std::runtime_error("Invalid configuration.  m_cfg.gamma_scale_corr_coeff_ecal should "
0049                              "have at least 3 parameters");
0050   }
0051   trace("gamma detection params:   max length={},   max width={},   max z={}", m_cfg.gammaMaxLength,
0052         m_cfg.gammaMaxWidth, m_gammaZMax);
0053 }
0054 /** calculates the correction for a given uncorrected total energy and a set of coefficients*/
0055 double FarForwardNeutralsReconstruction::calc_corr(double Etot, const std::vector<double>& coeffs) {
0056   return coeffs[0] + coeffs[1] / sqrt(Etot) + coeffs[2] / Etot;
0057 }
0058
0059 /**
0060      check that the cluster position is within the correct range,
0061      and that the sqrt(largest eigenvalue) is less than gamma_max_length,
0062      and that the sqrt(second largest eigenvalue) is less than gamma_max_width
0063   */
0064 bool FarForwardNeutralsReconstruction::isGamma(const edm4eic::Cluster& cluster) const {
0065   double l1 = sqrt(cluster.getShapeParameters(4)) * dd4hep::mm;
0066   double l2 = sqrt(cluster.getShapeParameters(5)) * dd4hep::mm;
0067   double l3 = sqrt(cluster.getShapeParameters(6)) * dd4hep::mm;
0068
0069   //z in the local coordinates
0070   double z = (cluster.getPosition().z * cos(m_cfg.globalToProtonRotation) +
0071               cluster.getPosition().x * sin(m_cfg.globalToProtonRotation)) *
0072              dd4hep::mm;
0073   trace("z recon = {}", z);
0074   trace("l1 = {}, l2 = {}, l3 = {}", l1, l2, l3);
0075   bool isZMoreThanMax = (z > m_gammaZMax);
0076   bool isLengthMoreThanMax =
0077       (l1 > m_cfg.gammaMaxLength || l2 > m_cfg.gammaMaxLength || l3 > m_cfg.gammaMaxLength);
0078   bool areWidthsMoreThanMax = static_cast<int>(l1 > m_cfg.gammaMaxWidth) +
0079                                   static_cast<int>(l2 > m_cfg.gammaMaxWidth) +
0080                                   static_cast<int>(l3 > m_cfg.gammaMaxWidth) >=
0081                               2;
0082   return !(isZMoreThanMax || isLengthMoreThanMax || areWidthsMoreThanMax);
0083 }
0084
0085 void FarForwardNeutralsReconstruction::process(
0086     const FarForwardNeutralsReconstruction::Input& input,
0087     const FarForwardNeutralsReconstruction::Output& output) const {
0088   const auto [clustersHcal] = input;
0089   auto [out_neutrals]       = output;
0090
0091   double Etot_hcal = 0;
0092   double Emax      = 0;
0093   edm4hep::Vector3f n_position;
0094   for (const auto& cluster : *clustersHcal) {
0095     double E = cluster.getEnergy();
0096
0097     if (isGamma(cluster)) {
0098       auto rec_part = out_neutrals->create();
0099       rec_part.setPDG(22);
0100       edm4hep::Vector3f position = cluster.getPosition();
0101       double corr                = calc_corr(E, m_cfg.gammaScaleCorrCoeffHcal);
0102       E                          = E / (1 + corr);
0103       double p                   = E;
0104       double r                   = edm4hep::utils::magnitude(position);
0105       edm4hep::Vector3f momentum = position * (p / r);
0106       rec_part.setEnergy(E);
0107       rec_part.setMomentum(momentum);
0108       rec_part.setReferencePoint(position);
0109       rec_part.setCharge(0);
0110       rec_part.setMass(0);
0111       rec_part.addToClusters(cluster);
0112       continue;
0113     }
0114
0115     Etot_hcal += E;
0116     if (E > Emax) {
0117       Emax       = E;
0118       n_position = cluster.getPosition();
0119     }
0120   }
0121
0122   double Etot                   = Etot_hcal;
0123   static const double m_neutron = m_particleSvc.particle(2112).mass;
0124   int n_neutrons                = 0;
0125   if (Etot > 0 && Emax > 0) {
0126     auto rec_part = out_neutrals->create();
0127     double corr   = calc_corr(Etot, m_cfg.neutronScaleCorrCoeffHcal);
0128     Etot_hcal     = Etot_hcal / (1 + corr);
0129     Etot          = Etot_hcal;
0130     rec_part.setEnergy(Etot);
0131     rec_part.setPDG(2112);
0132     double p                   = sqrt(Etot * Etot - m_neutron * m_neutron);
0133     double r                   = edm4hep::utils::magnitude(n_position);
0134     edm4hep::Vector3f momentum = n_position * (p / r);
0135     rec_part.setMomentum(momentum);
0136     rec_part.setReferencePoint(n_position);
0137     rec_part.setCharge(0);
0138     rec_part.setMass(m_neutron);
0139     for (const auto& cluster : *clustersHcal) {
0140       rec_part.addToClusters(cluster);
0141     }
0142     n_neutrons = 1;
0143   } else {
0144     n_neutrons = 0;
0145   }
0146   debug("Found {} neutron candidates", n_neutrons);
0147 }
0148 } // namespace eicrecon