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       if (m_cfg.neutronScaleCorrCoeffHcal.size() < 3) {
0030         error("Invalid configuration.  m_cfg.neutronScaleCorrCoeffHcal should have at least 3 parameters");
0031         throw std::runtime_error("Invalid configuration.  m_cfg.neutronScaleCorrCoeffHcal should have at least 3 parameters");
0032       }
0033       if (m_cfg.gammaScaleCorrCoeffHcal.size() < 3) {
0034         error("Invalid configuration.  m_cfg.gamma_scale_corr_coeff_ecal should have at least 3 parameters");
0035         throw std::runtime_error("Invalid configuration.  m_cfg.gamma_scale_corr_coeff_ecal should have at least 3 parameters");
0036       }
0037       trace("gamma detection params:   max length={},   max width={},   max z={}", m_cfg.gammaMaxLength, m_cfg.gammaMaxWidth,
0038             m_cfg.gammaZMax);
0039     }
0040     /** calculates the correction for a given uncorrected total energy and a set of coefficients*/
0041     double FarForwardNeutralsReconstruction::calc_corr(double Etot, const std::vector<double>& coeffs) const{
0042       return coeffs[0]+coeffs[1]/sqrt(Etot)+coeffs[2]/Etot;
0043     }
0044
0045   /**
0046      check that the cluster position is within the correct range,
0047      and that the sqrt(largest eigenvalue) is less than gamma_max_length,
0048      and that the sqrt(second largest eigenvalue) is less than gamma_max_width
0049   */
0050     bool FarForwardNeutralsReconstruction::isGamma(const edm4eic::Cluster& cluster) const{
0051       double l1=sqrt(cluster.getShapeParameters(4))*dd4hep::mm;
0052       double l2=sqrt(cluster.getShapeParameters(5))*dd4hep::mm;
0053       double l3=sqrt(cluster.getShapeParameters(6))*dd4hep::mm;
0054
0055       //z in the local coordinates
0056       double z= (cluster.getPosition().z*cos(m_cfg.globalToProtonRotation)+cluster.getPosition().x*sin(m_cfg.globalToProtonRotation))*dd4hep::mm;
0057       trace("z recon = {}", z);
0058       trace("l1 = {}, l2 = {}, l3 = {}", l1, l2, l3);
0059       bool isZMoreThanMax = (z > m_cfg.gammaZMax);
0060       bool isLengthMoreThanMax = (l1> m_cfg.gammaMaxLength || l2> m_cfg.gammaMaxLength || l3 > m_cfg.gammaMaxLength);
0061       bool areWidthsMoreThanMax = (l1> m_cfg.gammaMaxWidth)+(l2> m_cfg.gammaMaxWidth)+(l3 > m_cfg.gammaMaxWidth)>=2;
0062       return !(isZMoreThanMax || isLengthMoreThanMax || areWidthsMoreThanMax);
0063
0064     }
0065
0066
0067
0068     void FarForwardNeutralsReconstruction::process(const FarForwardNeutralsReconstruction::Input& input,
0069                       const FarForwardNeutralsReconstruction::Output& output) const {
0070       const auto [clustersHcal] = input;
0071       auto [out_neutrals] = output;
0072
0073       double Etot_hcal=0;
0074       double Emax=0;
0075       edm4hep::Vector3f n_position;
0076       for (const auto& cluster : *clustersHcal) {
0077           double E = cluster.getEnergy();
0078
0079           if(isGamma(cluster)){
0080             auto rec_part = out_neutrals->create();
0081             rec_part.setPDG(22);
0082             edm4hep::Vector3f position = cluster.getPosition();
0083             double corr=calc_corr(E,m_cfg.gammaScaleCorrCoeffHcal);
0084             E=E/(1+corr);
0085             double p = E;
0086             double r = edm4hep::utils::magnitude(position);
0087             edm4hep::Vector3f momentum = position * (p / r);
0088             rec_part.setEnergy(E);
0089             rec_part.setMomentum(momentum);
0090             rec_part.setReferencePoint(position);
0091             rec_part.setCharge(0);
0092             rec_part.setMass(0);
0093             rec_part.addToClusters(cluster);
0094             continue;
0095           }
0096
0097           Etot_hcal += E;
0098           if (E > Emax){
0099             Emax = E;
0100             n_position = cluster.getPosition();
0101           }
0102       }
0103
0104       double Etot=Etot_hcal;
0105       static const double m_neutron = m_particleSvc.particle(2112).mass;
0106       int n_neutrons;
0107       if (Etot > 0 && Emax > 0){
0108           auto rec_part = out_neutrals->create();
0109           double corr=calc_corr(Etot,m_cfg.neutronScaleCorrCoeffHcal);
0110           Etot_hcal=Etot_hcal/(1+corr);
0111           Etot=Etot_hcal;
0112           rec_part.setEnergy(Etot);
0113           rec_part.setPDG(2112);
0114           double p = sqrt(Etot*Etot-m_neutron*m_neutron);
0115           double r = edm4hep::utils::magnitude(n_position);
0116           edm4hep::Vector3f momentum = n_position * (p / r);
0117           rec_part.setMomentum(momentum);
0118           rec_part.setReferencePoint(n_position);
0119           rec_part.setCharge(0);
0120           rec_part.setMass(m_neutron);
0121           for (const auto& cluster : *clustersHcal){
0122             rec_part.addToClusters(cluster);
0123           }
0124           n_neutrons=1;
0125       } else {
0126         n_neutrons=0;
0127       }
0128       debug("Found {} neutron candidates", n_neutrons);
0129
0130     }
0131 }