Skip to content

Zdc pid #2740

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
sebouh137 wants to merge 6 commits into
base: main
Choose a base branch
from
Draft

Zdc pid #2740

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
e60bcad
first commit
sebouh137 Jun 9, 2026
1127c0c
Refactor far-forward neutral reconstruction for multiple calorimeters
baptistefraisse Jul 9, 2026
33ae708
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Jul 9, 2026
92f54d0
Merge branch 'main' into zdc_pid
sebouh137 Jul 13, 2026
44bfe02
git merge
sebouh137 Jul 13, 2026
d535c90
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Jul 13, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions compile.sh
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
if [[ -n $1 ]]; then
rm -rf build
fi
cmake -B build -S . -DCMAKE_INSTALL_PREFIX=install && cmake --build build -j10; cmake --install build
105 changes: 38 additions & 67 deletions src/algorithms/reco/FarForwardNeutralsReconstruction.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
#include <stdexcept>
#include <tuple>
#include <vector>
#include <string>

#include "FarForwardNeutralsReconstruction.h"

Expand Down Expand Up @@ -46,6 +47,8 @@ bool FarForwardNeutralsReconstruction::isGamma(const edm4eic::Cluster& cluster)
cluster.getPosition().x * sin(m_cfg.globalToProtonRotation)) *
dd4hep::mm;

double E = cluster.getEnergy();

trace("z recon = {}", z);
trace("l1 = {}, l2 = {}, l3 = {}", l1, l2, l3);

Expand All @@ -56,8 +59,10 @@ bool FarForwardNeutralsReconstruction::isGamma(const edm4eic::Cluster& cluster)
static_cast<int>(l2 > m_cfg.gammaMaxWidth) +
static_cast<int>(l3 > m_cfg.gammaMaxWidth) >=
2;

return !(isZMoreThanMax || isLengthMoreThanMax || areWidthsMoreThanMax);
// max number of hits for a photon shower for a given energy E is of the form a*E+b*sqrt(E)
bool hasMoreHitsThanMax =
cluster.getNhits() > m_cfg.gammaMaxNhitsCoeffLin * E + m_cfg.gammaMaxNhitsCoeffSqrt * sqrt(E);
return !(isZMoreThanMax || isLengthMoreThanMax || areWidthsMoreThanMax || hasMoreHitsThanMax);
}

double FarForwardNeutralsReconstruction::corrPower(double E, const std::vector<double>& coeffs) {
Expand All @@ -75,10 +80,8 @@ int FarForwardNeutralsReconstruction::processNeutralCalo(
edm4eic::ReconstructedParticleCollection* out_neutrals,
const std::vector<double>& gammaScaleCoeff, const std::vector<double>& neutronScaleCoeff,
bool canDetectGammas, bool canDetectNeutrons, const CorrFunc& gammaCorr,
const CorrFunc& neutronCorr, GammaMode gammaMode, double gammaLeaderFracMin, double clusterEmin,
NeutronMode neutronMode, bool associateAllClustersToNeutron) const {

(void)gammaLeaderFracMin;
const CorrFunc& neutronCorr, GammaMode gammaMode, double clusterEmin, NeutronMode neutronMode,
bool associateAllClustersToNeutron) const {

const double m_neutron = m_particleSvc.particle(2112).mass;

Expand Down Expand Up @@ -248,67 +251,35 @@ void FarForwardNeutralsReconstruction::process(
const FarForwardNeutralsReconstruction::Output& output) const {

// Unpacking
const auto [clustersHcal, clustersB0, clustersEcalEndcapP, clustersLFHCAL] = input;
auto [out_neutralsHcal, out_neutralsB0, out_neutralsEcalEndcapP, out_neutralsLFHCAL] = output;

// Global
int n_neutrons = 0;

// ZDC-Hcal
n_neutrons += processNeutralCalo(clustersHcal, out_neutralsHcal,
/*gammaScaleCoeff=*/m_cfg.gammaScaleCorrCoeffHcalZDC,
/*neutronScaleCoeff=*/m_cfg.neutronScaleCorrCoeffHcalZDC,
/*canDetectGammas=*/true,
/*canDetectNeutrons=*/true,
/*gammaCorr=*/corrPower,
/*neutronCorr=*/corrPower,
/*gammaMode=*/GammaMode::AllPassing,
/*gammaLeaderFracMin=*/0.0,
/*clusterEmin=*/m_cfg.clusterEminHcalZDC,
/*neutronMode=*/NeutronMode::SumAll,
/*associateAllClustersToNeutron=*/true);

// B0-Ecal
n_neutrons += processNeutralCalo(clustersB0, out_neutralsB0,
/*gammaScaleCoeff=*/m_cfg.gammaScaleCorrCoeffB0Ecal,
/*neutronScaleCoeff=*/m_cfg.neutronScaleCorrCoeffB0Ecal,
/*canDetectGammas=*/true,
/*canDetectNeutrons=*/false,
/*gammaCorr=*/corrPower,
/*neutronCorr=*/corrPower,
/*gammaMode=*/GammaMode::LeaderOnly,
/*gammaLeaderFracMin=*/0.0,
/*clusterEmin=*/m_cfg.clusterEminB0Ecal,
/*neutronMode=*/NeutronMode::None,
/*associateAllClustersToNeutron=*/false);

// EndcapP-Ecal
n_neutrons += processNeutralCalo(clustersEcalEndcapP, out_neutralsEcalEndcapP,
/*gammaScaleCoeff=*/m_cfg.gammaScaleCorrCoeffEcalEndcapP,
/*neutronScaleCoeff=*/m_cfg.neutronScaleCorrCoeffEcalEndcapP,
/*canDetectGammas=*/true,
/*canDetectNeutrons=*/false,
/*gammaCorr=*/corrPower,
/*neutronCorr=*/corrPower,
/*gammaMode=*/GammaMode::LeaderOnly,
/*gammaLeaderFracMin=*/0.0,
/*clusterEmin=*/m_cfg.clusterEminEcalEndcapP,
/*neutronMode=*/NeutronMode::None,
/*associateAllClustersToNeutron=*/false);

// LFHCAL
n_neutrons += processNeutralCalo(clustersLFHCAL, out_neutralsLFHCAL,
/*gammaScaleCoeff=*/m_cfg.gammaScaleCorrCoeffLFHCAL,
/*neutronScaleCoeff=*/m_cfg.neutronScaleCorrCoeffLFHCAL,
/*canDetectGammas=*/false,
/*canDetectNeutrons=*/true,
/*gammaCorr=*/corrPower,
/*neutronCorr=*/corrPower,
/*gammaMode=*/GammaMode::None,
/*gammaLeaderFracMin=*/0.0,
/*clusterEmin=*/m_cfg.clusterEminLFHCAL,
/*neutronMode=*/NeutronMode::LeaderOnly,
/*associateAllClustersToNeutron=*/false);
const auto [clusters] = input;
auto [out_neutrals] = output;

GammaMode gammaMode = GammaMode::None;
if (m_cfg.gammaMode == "LeaderOnly") {
gammaMode = GammaMode::LeaderOnly;
} else if (m_cfg.gammaMode == "AllPassing") {
gammaMode = GammaMode::AllPassing;
}

NeutronMode neutronMode = NeutronMode::None;
if (m_cfg.neutronMode == "SumAll") {
neutronMode = NeutronMode::SumAll;
} else if (m_cfg.neutronMode == "LeaderOnly") {
neutronMode = NeutronMode::LeaderOnly;
}

const int n_neutrons =
processNeutralCalo(clusters, out_neutrals,
/*gammaScaleCoeff=*/m_cfg.gammaScaleCorrCoeff,
/*neutronScaleCoeff=*/m_cfg.neutronScaleCorrCoeff,
/*canDetectGammas=*/m_cfg.canDetectGammas,
/*canDetectNeutrons=*/m_cfg.canDetectNeutrons,
/*gammaCorr=*/corrPower,
/*neutronCorr=*/corrPower,
/*gammaMode=*/gammaMode,
/*clusterEmin=*/m_cfg.clusterEmin,
/*neutronMode=*/neutronMode,
/*associateAllClustersToNeutron=*/m_cfg.associateAllClustersToNeutron);

debug("Found {} neutron candidates", n_neutrons);
}
Expand Down
27 changes: 6 additions & 21 deletions src/algorithms/reco/FarForwardNeutralsReconstruction.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,18 +21,8 @@
namespace eicrecon {

using FarForwardNeutralsReconstructionAlgorithm =

algorithms::Algorithm<
algorithms::Input<const edm4eic::ClusterCollection, // clusters ZDC-Hcal
const edm4eic::ClusterCollection, // clusters B0-Ecal
const edm4eic::ClusterCollection, // clusters EndcapP-Ecal
const edm4eic::ClusterCollection>, // clusters LFHCAL

algorithms::Output<
edm4eic::ReconstructedParticleCollection, // neutrons/gamma in ZDC-Hcal
edm4eic::ReconstructedParticleCollection, // neutrons/gamma in B0-Ecal
edm4eic::ReconstructedParticleCollection, // neutrons/gamma in EndcapP-Ecal
edm4eic::ReconstructedParticleCollection>>; // neutrons/gamma in LFHCAL
algorithms::Algorithm<algorithms::Input<const edm4eic::ClusterCollection>,
algorithms::Output<edm4eic::ReconstructedParticleCollection>>;
/**
* Reconstructs far-forward neutral candidates from multiple calorimeter cluster collections.
*
Expand All @@ -59,12 +49,8 @@ class FarForwardNeutralsReconstruction
FarForwardNeutralsReconstruction(std::string_view name)
: FarForwardNeutralsReconstructionAlgorithm{
name,

{"clustersHcal", "clustersB0", "clustersEcalEndCapP", "clustersLFHCAL"},

{"outputNeutralsHcal", "outputNeutralsB0", "outputNeutralsEcalEndCapP",
"outputNeutralsLFHCAL"},

{"clusters"},
{"outputNeutrals"},
"Convert EMCal and HCal clusters into neutron or photon candidates"} {}

void init() final;
Expand All @@ -89,9 +75,8 @@ class FarForwardNeutralsReconstruction
const std::vector<double>& gammaScaleCoeff,
const std::vector<double>& neutronScaleCoeff, bool canDetectGammas,
bool canDetectNeutrons, const CorrFunc& gammaCorr,
const CorrFunc& neutronCorr, GammaMode gammaMode,
double gammaLeaderFracMin, double clusterEmin, NeutronMode neutronMode,
bool associateAllClustersToNeutron) const;
const CorrFunc& neutronCorr, GammaMode gammaMode, double clusterEmin,
NeutronMode neutronMode, bool associateAllClustersToNeutron) const;
};

} // namespace eicrecon
55 changes: 29 additions & 26 deletions src/algorithms/reco/FarForwardNeutralsReconstructionConfig.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,39 +3,42 @@
#pragma once
#include <float.h>
#include <DD4hep/Detector.h>
#include <string>
#include <vector>

namespace eicrecon {

struct FarForwardNeutralsReconstructionConfig {
/** detector constant describing distance to the ZDC */
/** detector constant describing distance reference position */
std::string offsetPositionName = "HcalFarForwardZDC_SiPMonTile_r_pos";
/** Correction factors for neutrons in the Hcal (ZDC) */
std::vector<double> neutronScaleCorrCoeffHcalZDC = {2.4, 0.89};
/** Correction factors for gammas in the Hcal (ZDC) */
std::vector<double> gammaScaleCorrCoeffHcalZDC = {1.1, 0.98};
/** Correction factors for neutrons in the LFHCAL */
std::vector<double> neutronScaleCorrCoeffLFHCAL = {2.55, 0.95};
/** Correction factors for gammas in the LFHCAL */
std::vector<double> gammaScaleCorrCoeffLFHCAL = {0., 0.};
/** Correction factors for neutrons in the B0-Ecal */
std::vector<double> neutronScaleCorrCoeffB0Ecal = {0., 0.};
/** Correction factors for gammas in the B0-Ecal */
std::vector<double> gammaScaleCorrCoeffB0Ecal = {0.99, 1.14};
/** Correction factors for neutrons in the Endcap-Ecal */
std::vector<double> neutronScaleCorrCoeffEcalEndcapP = {0., 0.};
/** Correction factors for gammas in the Endcap-Ecal */
std::vector<double> gammaScaleCorrCoeffEcalEndcapP = {1.05, 1.01};
/** Cluster thresholds */
double clusterEminHcalZDC = 0.0; // GeV
double clusterEminB0Ecal = 1.0; // GeV
double clusterEminEcalEndcapP = 1.0; // GeV
double clusterEminLFHCAL = 7.0; // GeV

/** Correction factors */
std::vector<double> neutronScaleCorrCoeff = {0.0, 0.0};
std::vector<double> gammaScaleCorrCoeff = {0.0, 0.0};

/** Detector capabilities */
bool canDetectGammas = true;
bool canDetectNeutrons = true;

/** Reconstruction modes */
std::string gammaMode = "None";
std::string neutronMode = "None";

/** Cluster threshold */
double clusterEmin = 0.0; // GeV

/** Whether all non-gamma clusters are summed into one neutron candidate */
bool associateAllClustersToNeutron = false;

/** rotation from global to local coordinates */
double globalToProtonRotation = -0.025;
/** Neutron-photon separation in HcalFarForwardZDC */
double gammaZMaxOffset = 400;
double gammaMaxLength = 100;
double gammaMaxWidth = 12;

/** Neutron-photon separation used for ZDC */
double gammaZMaxOffset = 400;
double gammaMaxLength = 100;
double gammaMaxWidth = 12;
double gammaMaxNhitsCoeffLin = 0.3;
double gammaMaxNhitsCoeffSqrt = 30;
};

} // namespace eicrecon
41 changes: 21 additions & 20 deletions src/algorithms/reco/LambdaReconstruction.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -373,10 +373,6 @@ void LambdaReconstruction::process(const LambdaReconstruction::Input& input,
const double mL = std::sqrt(m2L);
const double dL = mL - m_lambda;

if (std::abs(dL) > m_cfg.lambdaMassWindow * m_lambda) {
continue;
}

const double pi0_term = pp.dmpi0_cand / (m_cfg.pi0Window * m_pi0);
const double lambda_term = dL / (m_cfg.lambdaMassWindow * m_lambda);

Expand Down Expand Up @@ -417,36 +413,41 @@ void LambdaReconstruction::process(const LambdaReconstruction::Input& input,
// --------------------------------------------------------------------------

auto better = [&](const LambdaCandidate& a, const LambdaCandidate& b) -> bool {
// Primary criterion: physical compatibility.
if (a.chi2 != b.chi2) {
return a.chi2 < b.chi2;
}
// Tie-breakers only: prefer ZDC information when candidates are similarly good.
if (a.n_cat != b.n_cat) {
return static_cast<int>(a.n_cat) < static_cast<int>(b.n_cat);
return a.n_cat == NeutronCategory::ZDC;
}
if (a.g_cat != b.g_cat) {
return static_cast<int>(a.g_cat) < static_cast<int>(b.g_cat);
}
if (a.chi2 != b.chi2) {
return a.chi2 < b.chi2;
}
if (a.pz != b.pz) {
return a.pz > b.pz;
}
return a.E > b.E;
};

int best_k = 0;
for (int k = 1; k < static_cast<int>(cands.size()); ++k) {
if (better(cands[k], cands[best_k])) {
best_k = k;
}
}
std::sort(cands.begin(), cands.end(), better);

@github-actions github-actions Bot Jul 13, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ modernize-use-ranges ⚠️
use a ranges version of this algorithm

Suggested change
std::sort(cands.begin(), cands.end(), better);
std::ranges::sort(cands,, better);


const int max_trials = cands.size();

const auto& best = cands[best_k];
for (int k = 0; k < max_trials; ++k) {
const auto& best = cands[k];

const auto& g1 = gamma_pool[best.g_i];
const auto& g2 = gamma_pool[best.g_j];
const auto& n =
(best.n_cat == NeutronCategory::ZDC) ? neutrons_zdc[best.n_idx] : neutrons_other[best.n_idx];
const auto& g1 = gamma_pool[best.g_i];
const auto& g2 = gamma_pool[best.g_j];
const auto& n = (best.n_cat == NeutronCategory::ZDC) ? neutrons_zdc[best.n_idx]
: neutrons_other[best.n_idx];

reconstruct_from_triplet(n, g1, g2, out_lambdas, out_decay_products);
if (!reconstruct_from_triplet(n, g1, g2, out_lambdas, out_decay_products)) {
continue;
}

break;
}
}

} // namespace eicrecon
Loading
Loading