algorithms/reco/FarForwardNeutronReconstruction.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2024 Sebouh Paul
0003
0004 #include <edm4eic/ClusterCollection.h>
0005 #include <edm4eic/ReconstructedParticleCollection.h>
0006 #include <edm4hep/Vector3f.h>
0007 #include <edm4hep/utils/vector_utils.h>
0008 #include <math.h>
0009 #include <gsl/pointers>
0010 #include <stdexcept>
0011 #include <vector>
0012
0013 #include "FarForwardNeutronReconstruction.h"
0014
0015 /**
0016  Creates a "neutron candidate" Reconstructed Particle consisting of all clusters from a
0017  specified hadronic calorimeter and (optionally) from a specified Electromagnetic calorimeter.
0018  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 factor is 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.  This form was chosen
0022  empirically based on the discrepancies in single-neutron MC simulations between the uncorrected
0023  reconstructed energies and the truth energies of the neutrons.  Separate correction factors
0024  are included for the Hcal and Ecal.
0025  */
0026
0027 namespace eicrecon {
0028
0029     void FarForwardNeutronReconstruction::init() {
0030       if (m_cfg.scale_corr_coeff_hcal.size() < 3) {
0031         error("Invalid configuration.  m_cfg.scale_corr_coeff_hcal should have at least 3 parameters");
0032         throw std::runtime_error("Invalid configuration.  m_cfg.scale_corr_coeff_hcal should have at least 3 parameters");
0033       }
0034       if (m_cfg.scale_corr_coeff_ecal.size() < 3) {
0035         error("Invalid configuration.  m_cfg.scale_corr_coeff_ecal should have at least 3 parameters");
0036         throw std::runtime_error("Invalid configuration.  m_cfg.scale_corr_coeff_ecal should have at least 3 parameters");
0037       }
0038     }
0039     /** calculates the correction for a given uncorrected total energy and a set of coefficients*/
0040     double FarForwardNeutronReconstruction::calc_corr(double Etot, const std::vector<double>& coeffs) const{
0041       return coeffs[0]+coeffs[1]/sqrt(Etot)+coeffs[2]/Etot;
0042     }
0043     void FarForwardNeutronReconstruction::process(const FarForwardNeutronReconstruction::Input& input,
0044                       const FarForwardNeutronReconstruction::Output& output) const {
0045       const auto [clustersHcal,clustersEcal] = input;
0046       auto [out_neutrons] = output;
0047
0048       double Etot_hcal=0, Etot_ecal=0;
0049       double Emax=0;
0050       edm4hep::Vector3f position;
0051       for (const auto& cluster : *clustersHcal) {
0052           double E = cluster.getEnergy();
0053           Etot_hcal += E;
0054           if (E > Emax){
0055             Emax = E;
0056             position = cluster.getPosition();
0057           }
0058       }
0059       for (const auto& cluster : *clustersEcal) {
0060           double E = cluster.getEnergy();
0061           Etot_ecal+=E;
0062       }
0063       double Etot=Etot_hcal+Etot_ecal;
0064       if (Etot > 0 && Emax > 0){
0065           auto rec_part = out_neutrons->create();
0066           double corr=calc_corr(Etot,m_cfg.scale_corr_coeff_hcal);
0067           Etot_hcal=Etot_hcal/(1+corr);
0068           corr=calc_corr(Etot,m_cfg.scale_corr_coeff_ecal);
0069           Etot_ecal=Etot_ecal/(1+corr);
0070           Etot=Etot_hcal+Etot_ecal;
0071           rec_part.setEnergy(Etot);
0072           rec_part.setPDG(2112);
0073           double p = sqrt(Etot*Etot-m_neutron*m_neutron);
0074           double r = edm4hep::utils::magnitude(position);
0075           edm4hep::Vector3f momentum = position * (p / r);
0076           rec_part.setMomentum(momentum);
0077           rec_part.setCharge(0);
0078           rec_part.setMass(m_neutron);
0079           for (const auto& cluster : *clustersHcal){
0080             rec_part.addToClusters(cluster);
0081           }
0082           for (const auto& cluster : *clustersEcal){
0083             rec_part.addToClusters(cluster);
0084           }
0085       }
0086         //m_log->debug("Found {} neutron candidates", out_neutrons->size());
0087
0088     }
0089 }