Lambda reconstruction in ZDC

src/algorithms/calorimetry/HEXPLIT.cc

@@ -82,6 +82,7 @@ void HEXPLIT::process(const HEXPLIT::Input& input,
   auto [subcellHits] = output;
 
   double MIP=m_cfg.MIP/dd4hep::GeV;
+  double delta=m_cfg.delta_in_MIPs*MIP;
   double Emin=m_cfg.Emin_in_MIPs*MIP;
   double tmax=m_cfg.tmax/dd4hep::ns;
 
@@ -125,7 +126,7 @@ void HEXPLIT::process(const HEXPLIT::Input& input,
     }
     double weights[SUBCELLS];
     for(int k=0; k<NEIGHBORS; k++){
-      Eneighbors[k]=std::max(Eneighbors[k],MIP);
+      Eneighbors[k]=std::max(Eneighbors[k],delta);
     }
     double sum_weights=0;
     for(int k=0; k<SUBCELLS; k++){

src/algorithms/calorimetry/HEXPLITConfig.h

@@ -8,6 +8,7 @@ namespace eicrecon {
   struct HEXPLITConfig {
     double                   MIP{472.*dd4hep::keV};
     double                   Emin_in_MIPs{0.1};
+    double                   delta_in_MIPs{0.01};
     double                   tmax{325*dd4hep::ns};
   };
 

src/algorithms/pid_lut/PIDLookup.cc

@@ -5,6 +5,7 @@
 #include <edm4eic/MCRecoParticleAssociationCollection.h>
 #include <edm4eic/ReconstructedParticleCollection.h>
 #include <edm4hep/MCParticleCollection.h>
+#include <edm4hep/Vector3f.h>
 #include <edm4hep/utils/vector_utils.h>
 #include <fmt/core.h>
 #include <podio/podioVersion.h>

src/algorithms/reco/FarForwardLambdaReconstruction.cc

@@ -0,0 +1,157 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Sebouh Paul
+
+#include <Evaluator/DD4hepUnits.h>
+#include <TVector3.h>
+#include <edm4eic/ReconstructedParticleCollection.h>
+#include <edm4hep/Vector3f.h>
+#include <edm4hep/utils/vector_utils.h>
+#include <fmt/core.h>
+#include <math.h>
+#include <cstddef>
+#include <gsl/pointers>
+#include <vector>
+
+#include "FarForwardLambdaReconstruction.h"
+#include "TLorentzVector.h"
+/**
+Creates Lambda candidates from a neutron and two photons from a pi0 decay
+ */
+
+namespace eicrecon {
+
+    void FarForwardLambdaReconstruction::init() {
+
+    }
+
+  /* converts one type of vector format to another */
+  void toTVector3(TVector3& v1, const edm4hep::Vector3f& v2){
+    v1.SetXYZ(v2.x,v2.y,v2.z);
+  }
+
+    void FarForwardLambdaReconstruction::process(const FarForwardLambdaReconstruction::Input& input,
+                      const FarForwardLambdaReconstruction::Output& output) const {
+      const auto [neutrals] = input;
+      auto [out_lambdas] = output;
+      std::vector<edm4eic::ReconstructedParticle> neutrons={};
+      std::vector<edm4eic::ReconstructedParticle> gammas={};
+      for (auto part: *neutrals){
+        if (part.getPDG()==2112){
+          neutrons.push_back(part);
+        }
+        if (part.getPDG()==22){
+          gammas.push_back(part);
+        }
+      }
+
+      if (neutrons.size()<1 || gammas.size()<2)
+        return;
+
+      static const double m_neutron = m_particleSvc.particle(2112).mass;
+      static const double m_pi0 = m_particleSvc.particle(111).mass;
+      static const double m_lambda = m_particleSvc.particle(3122).mass;
+
+      for (std::size_t i_n=0; i_n<neutrons.size(); i_n++){
+        for (std::size_t i_1=0; i_1<gammas.size()-1; i_1++){
+          for (std::size_t i_2=i_1+1; i_2<gammas.size(); i_2++){
+
+            double En=neutrons[i_n].getEnergy();
+            double pn=sqrt(En*En-m_neutron*m_neutron);
+            double E1=gammas[i_1].getEnergy();
+            double E2=gammas[i_2].getEnergy();
+            TVector3 xn;
+            TVector3 x1;
+            TVector3 x2;
+
+            toTVector3(xn,neutrons[0].getReferencePoint()*dd4hep::mm);
+            toTVector3(x1,gammas[0].getReferencePoint()*dd4hep::mm);
+            toTVector3(x2,gammas[1].getReferencePoint()*dd4hep::mm);
+
+            xn.RotateY(-m_cfg.globalToProtonRotation);
+            x1.RotateY(-m_cfg.globalToProtonRotation);
+            x2.RotateY(-m_cfg.globalToProtonRotation);
+
+            debug("nx recon = {}, g1x recon = {}, g2x recon = {}", xn.X(), x1.X(), x2.X());
+            debug("nz recon = {}, g1z recon = {}, g2z recon = {}, z face = {}", xn.Z(), x1.Z(), x2.Z(), m_cfg.zMax);
+
+            TVector3 vtx(0,0,0);
+            double f=0;
+            double df=0.5;
+            double theta_open_expected=2*asin(m_pi0/(2*sqrt(E1*E2)));
+            TLorentzVector n, g1, g2, lambda;
+            for(int i = 0; i<m_cfg.iterations; i++){
+              n={pn*(xn-vtx).Unit(), En};
+              g1={E1*(x1-vtx).Unit(), E1};
+              g2={E2*(x2-vtx).Unit(), E2};
+              double theta_open=g1.Angle(g2.Vect());
+              lambda=n+g1+g2;
+              if (theta_open>theta_open_expected) {
+                f-=df;
+              } else if (theta_open<theta_open_expected) {
+                f+=df;
+              }
+
+              vtx=lambda.Vect()*(f*m_cfg.zMax/lambda.Z());
+              df/=2;
+            }
+
+            double mass_rec=lambda.M();
+            if (abs(mass_rec-m_lambda)>m_cfg.lambdaMaxMassDev)
+              continue;
+
+            // rotate everything back to the lab coordinates.
+            vtx.RotateY(m_cfg.globalToProtonRotation);
+            lambda.RotateY(m_cfg.globalToProtonRotation);
+            n.RotateY(m_cfg.globalToProtonRotation);
+            g1.RotateY(m_cfg.globalToProtonRotation);
+            g2.RotateY(m_cfg.globalToProtonRotation);
+
+            auto b=-lambda.BoostVector();
+            n.Boost(b);
+            g1.Boost(b);
+            g2.Boost(b);
+
+            //convert vertex to mm:
+            vtx=vtx*(1/dd4hep::mm);
+
+            auto rec_part = out_lambdas->create();
+            rec_part.setPDG(3122);
+
+            rec_part.setEnergy(lambda.E());
+            rec_part.setMomentum({static_cast<float>(lambda.X()), static_cast<float>(lambda.Y()), static_cast<float>(lambda.Z())});
+            rec_part.setReferencePoint({static_cast<float>(vtx.X()), static_cast<float>(vtx.Y()), static_cast<float>(vtx.Z())});
+            rec_part.setCharge(0);
+            rec_part.setMass(mass_rec);
+
+            auto neutron = out_lambdas->create();
+            neutron.setPDG(2112);
+            neutron.setEnergy(n.E());
+            neutron.setMomentum({static_cast<float>(n.X()), static_cast<float>(n.Y()), static_cast<float>(n.Z())});
+            neutron.setReferencePoint({static_cast<float>(vtx.X()), static_cast<float>(vtx.Y()), static_cast<float>(vtx.Z())});
+            neutron.setCharge(0);
+            neutron.setMass(m_neutron);
+            rec_part.addToParticles(neutron);
+
+            auto gamma1 = out_lambdas->create();
+            gamma1.setPDG(22);
+            gamma1.setEnergy(g1.E());
+            gamma1.setMomentum({static_cast<float>(g1.X()), static_cast<float>(g1.Y()), static_cast<float>(g1.Z())});
+            gamma1.setReferencePoint({static_cast<float>(vtx.X()), static_cast<float>(vtx.Y()), static_cast<float>(vtx.Z())});
+            gamma1.setCharge(0);
+            gamma1.setMass(0);
+            rec_part.addToParticles(gamma1);
+
+            auto gamma2 = out_lambdas->create();
+            gamma2.setPDG(22);
+            gamma2.setEnergy(g2.E());
+            gamma2.setMomentum({static_cast<float>(g2.X()), static_cast<float>(g2.Y()), static_cast<float>(g2.Z())});
+            gamma2.setReferencePoint({static_cast<float>(vtx.X()), static_cast<float>(vtx.Y()), static_cast<float>(vtx.Z())});
+            gamma2.setCharge(0);
+            gamma2.setMass(0);
+            rec_part.addToParticles(gamma2);
+            continue;
+          }
+        }
+      }
+    }
+}

src/algorithms/reco/FarForwardLambdaReconstruction.h

@@ -0,0 +1,45 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Sebouh Paul
+
+#pragma once
+#include <algorithms/algorithm.h>
+#include <edm4eic/ReconstructedParticleCollection.h>
+#include <spdlog/logger.h>
+#include <memory>
+#include <string>                                 // for basic_string
+#include <string_view>                            // for string_view
+
+#include "algorithms/interfaces/ParticleSvc.h"
+#include "algorithms/interfaces/WithPodConfig.h"
+#include "algorithms/reco/FarForwardLambdaReconstructionConfig.h"
+
+namespace eicrecon {
+
+using FarForwardLambdaReconstructionAlgorithm = algorithms::Algorithm<
+   algorithms::Input<
+      const edm4eic::ReconstructedParticleCollection
+    >,
+  /*output collection contains the lambda candidates and their decay products in the CM frame*/
+    algorithms::Output<
+      edm4eic::ReconstructedParticleCollection
+    >
+    >;
+    class FarForwardLambdaReconstruction :
+       public FarForwardLambdaReconstructionAlgorithm,
+       public WithPodConfig<FarForwardLambdaReconstructionConfig> {
+       public:
+         FarForwardLambdaReconstruction(std::string_view name)
+                  : FarForwardLambdaReconstructionAlgorithm{name,
+                                                            {"inputNeutrals"},
+                                                            {"outputLambdaAndDecayProducts"},
+                                        "Reconstructs lambda candidates and their decay products from the reconstructed neutrons and photons"} {}
+
+         void init() final;
+         void process(const Input&, const Output&) const final;
+
+    private:
+        std::shared_ptr<spdlog::logger> m_log;
+       const algorithms::ParticleSvc& m_particleSvc = algorithms::ParticleSvc::instance();
+
+    };
+} // namespace eicrecon

src/algorithms/reco/FarForwardLambdaReconstructionConfig.h

@@ -0,0 +1,21 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Sebouh Paul
+#pragma once
+#include <float.h>
+#include <DD4hep/Detector.h>
+
+
+namespace eicrecon {
+
+  struct FarForwardLambdaReconstructionConfig {
+    /** transformation from global coordinates to proton-frame coordinates*/
+    double globalToProtonRotation=-0.025;
+    /** distance to the ZDC */
+    double zMax=35800*dd4hep::mm;
+    /** maximum deviation between reconstructed mass and PDG mass */
+    double lambdaMaxMassDev=0.030*dd4hep::GeV;
+    /** number of iterations for the IDOLA algorithm */
+    int iterations=10;
+  };
+
+} // eicrecon

src/algorithms/reco/FarForwardNeutralsReconstruction.cc

@@ -0,0 +1,131 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Sebouh Paul
+#include <Evaluator/DD4hepUnits.h>
+#include <edm4eic/ClusterCollection.h>
+#include <edm4eic/ReconstructedParticleCollection.h>
+#include <edm4hep/Vector3f.h>
+#include <edm4hep/utils/vector_utils.h>
+#include <fmt/core.h>
+#include <cmath>
+#include <gsl/pointers>
+#include <stdexcept>
+#include <vector>
+
+#include "FarForwardNeutralsReconstruction.h"
+
+/**
+ 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.
+ 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
+ times a correction factor, and its direction is the direction from the origin to the position
+ 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),
+ 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
+ empirically based on the discrepancies in single-neutron and single-photon MC simulations between the uncorrected
+ reconstructed energies and the truth energies of the neutrons.  The parameter values should be co-tuned with those of the clustering algorithm being used.
+ */
+
+namespace eicrecon {
+
+    void FarForwardNeutralsReconstruction::init() {
+      if (m_cfg.neutronScaleCorrCoeffHcal.size() < 3) {
+        error("Invalid configuration.  m_cfg.neutronScaleCorrCoeffHcal should have at least 3 parameters");
+        throw std::runtime_error("Invalid configuration.  m_cfg.neutronScaleCorrCoeffHcal should have at least 3 parameters");
+      }
+      if (m_cfg.gammaScaleCorrCoeffHcal.size() < 3) {
+        error("Invalid configuration.  m_cfg.gamma_scale_corr_coeff_ecal should have at least 3 parameters");
+        throw std::runtime_error("Invalid configuration.  m_cfg.gamma_scale_corr_coeff_ecal should have at least 3 parameters");
+      }
+      trace("gamma detection params:   max length={},   max width={},   max z={}", m_cfg.gammaMaxLength, m_cfg.gammaMaxWidth,
+            m_cfg.gammaZMax);
+    }
+    /** calculates the correction for a given uncorrected total energy and a set of coefficients*/
+    double FarForwardNeutralsReconstruction::calc_corr(double Etot, const std::vector<double>& coeffs) const{
+      return coeffs[0]+coeffs[1]/sqrt(Etot)+coeffs[2]/Etot;
+    }
+
+  /**
+     check that the cluster position is within the correct range,
+     and that the sqrt(largest eigenvalue) is less than gamma_max_length,
+     and that the sqrt(second largest eigenvalue) is less than gamma_max_width
+  */
+    bool FarForwardNeutralsReconstruction::isGamma(const edm4eic::Cluster& cluster) const{
+      double l1=sqrt(cluster.getShapeParameters(4))*dd4hep::mm;
+      double l2=sqrt(cluster.getShapeParameters(5))*dd4hep::mm;
+      double l3=sqrt(cluster.getShapeParameters(6))*dd4hep::mm;
+
+      //z in the local coordinates
+      double z= (cluster.getPosition().z*cos(m_cfg.globalToProtonRotation)+cluster.getPosition().x*sin(m_cfg.globalToProtonRotation))*dd4hep::mm;
+      trace("z recon = {}", z);
+      trace("l1 = {}, l2 = {}, l3 = {}", l1, l2, l3);
+      bool isZMoreThanMax = (z > m_cfg.gammaZMax);
+      bool isLengthMoreThanMax = (l1> m_cfg.gammaMaxLength || l2> m_cfg.gammaMaxLength || l3 > m_cfg.gammaMaxLength);
+      bool areWidthsMoreThanMax = (l1> m_cfg.gammaMaxWidth)+(l2> m_cfg.gammaMaxWidth)+(l3 > m_cfg.gammaMaxWidth)>=2;
+      return !(isZMoreThanMax || isLengthMoreThanMax || areWidthsMoreThanMax);
+
+    }
+
+
+
+    void FarForwardNeutralsReconstruction::process(const FarForwardNeutralsReconstruction::Input& input,
+                      const FarForwardNeutralsReconstruction::Output& output) const {
+      const auto [clustersHcal] = input;
+      auto [out_neutrals] = output;
+
+      double Etot_hcal=0;
+      double Emax=0;
+      edm4hep::Vector3f n_position;
+      for (const auto& cluster : *clustersHcal) {
+          double E = cluster.getEnergy();
+
+          if(isGamma(cluster)){
+            auto rec_part = out_neutrals->create();
+            rec_part.setPDG(22);
+            edm4hep::Vector3f position = cluster.getPosition();
+            double corr=calc_corr(E,m_cfg.gammaScaleCorrCoeffHcal);
+            E=E/(1+corr);
+            double p = E;
+            double r = edm4hep::utils::magnitude(position);
+            edm4hep::Vector3f momentum = position * (p / r);
+            rec_part.setEnergy(E);
+            rec_part.setMomentum(momentum);
+            rec_part.setReferencePoint(position);
+            rec_part.setCharge(0);
+            rec_part.setMass(0);
+            rec_part.addToClusters(cluster);
+            continue;
+          }
+
+          Etot_hcal += E;
+          if (E > Emax){
+            Emax = E;
+            n_position = cluster.getPosition();
+          }
+      }
+
+      double Etot=Etot_hcal;
+      static const double m_neutron = m_particleSvc.particle(2112).mass;
+      int n_neutrons;
+      if (Etot > 0 && Emax > 0){
+          auto rec_part = out_neutrals->create();
+          double corr=calc_corr(Etot,m_cfg.neutronScaleCorrCoeffHcal);
+          Etot_hcal=Etot_hcal/(1+corr);
+          Etot=Etot_hcal;
+          rec_part.setEnergy(Etot);
+          rec_part.setPDG(2112);
+          double p = sqrt(Etot*Etot-m_neutron*m_neutron);
+          double r = edm4hep::utils::magnitude(n_position);
+          edm4hep::Vector3f momentum = n_position * (p / r);
+          rec_part.setMomentum(momentum);
+          rec_part.setReferencePoint(n_position);
+          rec_part.setCharge(0);
+          rec_part.setMass(m_neutron);
+          for (const auto& cluster : *clustersHcal){
+            rec_part.addToClusters(cluster);
+          }
+          n_neutrons=1;
+      } else {
+        n_neutrons=0;
+      }
+      debug("Found {} neutron candidates", n_neutrons);
+
+    }
+}