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,139 @@
+// 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 <math.h>
+#include <gsl/pointers>
+
+
+#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 [neutrons, gammas] = input;
+      auto [out_lambdas, out_decay_products] = output;
+
+      if (neutrons->size()!=1 || gammas->size()!=2)
+        return;
+
+      double En=(*neutrons)[0].getEnergy();
+      double pn=sqrt(En*En-m_neutron*m_neutron);
+      double E1=(*gammas)[0].getEnergy();
+      double E2=(*gammas)[1].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.rot_y);
+      x1.RotateY(-m_cfg.rot_y);
+      x2.RotateY(-m_cfg.rot_y);
+
+
+      /*std::cout << "nx recon " << xn.X() << " g1x recon " << x1.X() << " g2x recon " << x2.X() << std::endl;
+
+        std::cout << "nz recon " << xn.Z() << " g1z recon " << x1.Z() << " g2z recon " << x2.Z() << " z face=" << m_cfg.zmax<< std::endl;*/
+
+
+      TVector3 vtx(0,0,0);
+      double f=0;
+      double df=0.5;
+      double theta_open_expected=2*asin(m_pi0/(2*sqrt(E1*E2)));
+      for(int i = 0; i<m_cfg.iterations; i++){
+        TLorentzVector n(pn*(xn-vtx).Unit(), En);
+        TLorentzVector g1(E1*(x1-vtx).Unit(), E1);
+        TLorentzVector g2(E2*(x2-vtx).Unit(), E2);
+        double theta_open=g1.Angle(g2.Vect());
+        TLorentzVector 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;
+
+        if (i==m_cfg.iterations-1){
+          double mass_rec=lambda.M();
+          if (abs(mass_rec-m_lambda)>m_cfg.lambda_max_mass_dev)
+            return;
+
+          // rotate everything back to the lab coordinates.
+          vtx.RotateY(m_cfg.rot_y);
+          lambda.RotateY(m_cfg.rot_y);
+          n.RotateY(m_cfg.rot_y);
+          g1.RotateY(m_cfg.rot_y);
+          g2.RotateY(m_cfg.rot_y);
+
+          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);
+          //edm4hep::Vector3f position(vtx.X(), vtx.Y(), vtx.Z());
+          //edm4hep::Vector3f momentum(lambda.X(), lambda.Y(), lambda.Z());
+
+          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);
+
+          rec_part = out_decay_products->create();
+          rec_part.setPDG(2112);
+          rec_part.setEnergy(n.E());
+          rec_part.setMomentum({static_cast<float>(n.X()), static_cast<float>(n.Y()), static_cast<float>(n.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(m_neutron);
+
+          rec_part = out_decay_products->create();
+          rec_part.setPDG(22);
+          rec_part.setEnergy(g1.E());
+          rec_part.setMomentum({static_cast<float>(g1.X()), static_cast<float>(g1.Y()), static_cast<float>(g1.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(0);
+
+          rec_part = out_decay_products->create();
+          rec_part.setPDG(22);
+          rec_part.setEnergy(g2.E());
+          rec_part.setMomentum({static_cast<float>(g2.X()), static_cast<float>(g2.Y()), static_cast<float>(g2.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(0);
+
+          return;
+
+        }
+      }
+    }
+}

src/algorithms/reco/FarForwardLambdaReconstruction.h

@@ -0,0 +1,51 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Sebouh Paul
+
+#pragma once
+#include <DD4hep/Detector.h>
+#include <algorithms/algorithm.h>
+#include <algorithms/geo.h>
+#include <edm4eic/ReconstructedParticleCollection.h>
+#include <spdlog/logger.h>
+#include <gsl/pointers>
+#include <memory>
+#include <string>                                 // for basic_string
+#include <string_view>                            // for string_view
+
+#include "algorithms/interfaces/WithPodConfig.h"
+#include "algorithms/reco/FarForwardLambdaReconstructionConfig.h"
+
+namespace eicrecon {
+
+using FarForwardLambdaReconstructionAlgorithm = algorithms::Algorithm<
+   algorithms::Input<
+      const edm4eic::ReconstructedParticleCollection,
+      const edm4eic::ReconstructedParticleCollection
+    >,
+    algorithms::Output<
+      edm4eic::ReconstructedParticleCollection,
+      edm4eic::ReconstructedParticleCollection
+    >
+    >;
+    class FarForwardLambdaReconstruction :
+       public FarForwardLambdaReconstructionAlgorithm,
+       public WithPodConfig<FarForwardLambdaReconstructionConfig> {
+       public:
+         FarForwardLambdaReconstruction(std::string_view name)
+                  : FarForwardLambdaReconstructionAlgorithm{name,
+                                                            {"inputNeutrons", "inputPhotons"},
+                                                            {"outputLambda", "outputLambdaDecayProducts"},
+                                        "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;
+        double m_neutron{0.93956542052};
+      double m_pi0{0.1349768277676847};
+      double m_lambda{1.115683138712051};
+        const dd4hep::Detector* m_detector{algorithms::GeoSvc::instance().detector()};
+
+    };
+} // 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 rot_y=-0.025;
+    /** distance to the ZDC */
+    double zmax=35800*dd4hep::mm;
+    /** maximum deviation between reconstructed mass and PDG mass */
+    double lambda_max_mass_dev=0.030*dd4hep::GeV;
+    /** number of iterations for the IDOLA algorithm */
+    int iterations=10;
+  };
+
+} // eicrecon

src/algorithms/reco/FarForwardNeutralsReconstruction.cc

@@ -0,0 +1,133 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Sebouh Paul
+
+#include <edm4eic/ClusterCollection.h>
+#include <edm4eic/ReconstructedParticleCollection.h>
+#include <edm4hep/Vector3f.h>
+#include <edm4hep/utils/vector_utils.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.n_scale_corr_coeff_hcal.size() < 3) {
+        error("Invalid configuration.  m_cfg.n_scale_corr_coeff_hcal should have at least 3 parameters");
+        throw std::runtime_error("Invalid configuration.  m_cfg.n_scale_corr_coeff_hcal should have at least 3 parameters");
+      }
+      if (m_cfg.gamma_scale_corr_coeff_hcal.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");
+      }
+    }
+    /** 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.rot_y)+cluster.getPosition().x*sin(m_cfg.rot_y))*dd4hep::mm;
+
+      /*std::cout << "zcut=" << m_cfg.gamma_zmax << " z recon=" << z << std::endl;
+      std::cout << "l1=" << l1 << " l2=" << l2 << " l3=" << l3 << std::endl;
+      std::cout << "max length=" << m_cfg.gamma_max_length << " max width=" << m_cfg.gamma_max_width << std::endl;*/
+      if (z> m_cfg.gamma_zmax)
+        return false;
+      if (l1> m_cfg.gamma_max_length || l2> m_cfg.gamma_max_length || l3 > m_cfg.gamma_max_length)
+        return false;
+      if ((l1> m_cfg.gamma_max_width) + (l2> m_cfg.gamma_max_width) + (l3 > m_cfg.gamma_max_width)>2)
+        return false;
+      return true;
+
+    }
+
+
+
+    void FarForwardNeutralsReconstruction::process(const FarForwardNeutralsReconstruction::Input& input,
+                      const FarForwardNeutralsReconstruction::Output& output) const {
+      const auto [clustersHcal] = input;
+      auto [out_neutrons, out_gammas] = 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_gammas->create();
+            rec_part.setPDG(22);
+            edm4hep::Vector3f position = cluster.getPosition();
+            double corr=calc_corr(E,m_cfg.gamma_scale_corr_coeff_hcal);
+            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;
+      if (Etot > 0 && Emax > 0){
+          auto rec_part = out_neutrons->create();
+          double corr=calc_corr(Etot,m_cfg.n_scale_corr_coeff_hcal);
+          Etot_hcal=Etot_hcal/(1+corr);
+          //corr=calc_corr(Etot,m_cfg.scale_corr_coeff_ecal);
+          //Etot_ecal=Etot_ecal/(1+corr);
+          Etot=Etot_hcal;//+Etot_ecal;
+          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);
+          }
+          /*for (const auto& cluster : *clustersEcal){
+            rec_part.addToClusters(cluster);
+            }*/
+      }
+        //m_log->debug("Found {} neutron candidates", out_neutrons->size());
+
+    }
+}