From aa9beab077e6b1912798a118a08584c1ac554cdd Mon Sep 17 00:00:00 2001
From: Wouter Deconinck <wdconinc@gmail.com>
Date: Wed, 7 Aug 2024 15:26:09 -0500
Subject: [PATCH] fix: avoid division by zero in neutron identification when no
 Hcal cluster (#1548)

### Briefly, what does this PR introduce?
When there is no Hcal cluster found, then the position is still
initialized to zero. Setting the momentum components based on a
normalized position fails with a division by zero, see e.g.
https://github.com/eic/EICrecon/actions/runs/10237639088/job/28321546895#step:8:6883.

This PR requires a cluster in the Hcal for there to be a neutron
identification with a momentum direction given by that cluster.

This PR also modifies some of the calls to use a vector type, which is
semantically cleaner as it reduces the number of independent variables.

### What kind of change does this PR introduce?
- [x] Bug fix (issue: division by zero)
- [ ] New feature (issue #__)
- [ ] Documentation update
- [ ] Other: __

### Please check if this PR fulfills the following:
- [ ] Tests for the changes have been added
- [ ] Documentation has been added / updated
- [x] Changes have been communicated to collaborators: @sebouh137
@ruse-traveler

### Does this PR introduce breaking changes? What changes might users
need to make to their code?
No.

### Does this PR change default behavior?
No.
---
 .../reco/FarForwardNeutronReconstruction.cc   | 27 ++++++++-----------
 1 file changed, 11 insertions(+), 16 deletions(-)

diff --git a/src/algorithms/reco/FarForwardNeutronReconstruction.cc b/src/algorithms/reco/FarForwardNeutronReconstruction.cc
index e22eff025..14b63a504 100644
--- a/src/algorithms/reco/FarForwardNeutronReconstruction.cc
+++ b/src/algorithms/reco/FarForwardNeutronReconstruction.cc
@@ -4,6 +4,7 @@
 #include <edm4eic/ClusterCollection.h>
 #include <edm4eic/ReconstructedParticleCollection.h>
 #include <edm4hep/Vector3f.h>
+#include <edm4hep/utils/vector_utils.h>
 #include <math.h>
 #include <gsl/pointers>
 #include <stdexcept>
@@ -46,17 +47,13 @@ namespace eicrecon {
 
       double Etot_hcal=0, Etot_ecal=0;
       double Emax=0;
-      double x=0;
-      double y=0;
-      double z=0;
+      edm4hep::Vector3f position;
       for (const auto& cluster : *clustersHcal) {
           double E = cluster.getEnergy();
-          Etot_hcal+=E;
-          if(E>Emax){
-            Emax=E;
-            x=cluster.getPosition().x;
-            y=cluster.getPosition().y;
-            z=cluster.getPosition().z;
+          Etot_hcal += E;
+          if (E > Emax){
+            Emax = E;
+            position = cluster.getPosition();
           }
       }
       for (const auto& cluster : *clustersEcal) {
@@ -64,7 +61,7 @@ namespace eicrecon {
           Etot_ecal+=E;
       }
       double Etot=Etot_hcal+Etot_ecal;
-      if (Etot>0){
+      if (Etot > 0 && Emax > 0){
           auto rec_part = out_neutrons->create();
           double corr=calc_corr(Etot,m_cfg.scale_corr_coeff_hcal);
           Etot_hcal=Etot_hcal/(1+corr);
@@ -73,12 +70,10 @@ namespace eicrecon {
           Etot=Etot_hcal+Etot_ecal;
           rec_part.setEnergy(Etot);
           rec_part.setPDG(2112);
-          double p=sqrt(Etot*Etot-m_neutron*m_neutron);
-          double r=sqrt(x*x+y*y+z*z);
-          double px=p*x/r;
-          double py=p*y/r;
-          double pz=p*z/r;
-          rec_part.setMomentum({(float)px, (float)py, (float)pz});
+          double p = sqrt(Etot*Etot-m_neutron*m_neutron);
+          double r = edm4hep::utils::magnitude(position);
+          edm4hep::Vector3f momentum = position * (p / r);
+          rec_part.setMomentum(momentum);
           rec_part.setCharge(0);
           rec_part.setMass(m_neutron);
           for (const auto& cluster : *clustersHcal){