HNL hadronic and electroweak decays by austinschneider · Pull Request #110 · Harvard-Neutrino/SIREN

austinschneider · 2026-04-28T04:29:10Z

Stack: PR 12 of 14

This PR is part of a 14-PR stack decomposing dev/HNL_DIS into reviewable chunks.

Base: feat/hnl-decays-2body
Head: feat/hnl-decays-hadronic-ew

Merge order: feat/hnl-decays-2body should land before this PR.

Squashed contents

Squashed diff for paths:

projects/interactions/private/HNLDecay.cxx
projects/interactions/public/SIREN/interactions/HNLDecay.h
projects/interactions/private/pybindings/HNLDecay.h
projects/interactions/private/ElectroweakDecay.cxx
projects/interactions/public/SIREN/interactions/ElectroweakDecay.h
projects/interactions/private/HNLDipoleDecay.cxx
projects/interactions/public/SIREN/interactions/HNLDipoleDecay.h
projects/interactions/private/pybindings/HNLDipoleDecay.h
projects/interactions/private/pybindings/interactions.cxx

Source commits on dev/HNL_DIS_clean that contributed:
197bf14 (Nicholas Kamp): fixing build errors from the merge
1bbd678 (Nicholas Kamp): fix memory leak
2219939 (nickkamp1): first step at EW decay
285285f (Nicholas Kamp): kinematics fix to 2body decays
3625e81 (Nicholas Kamp): fix build errors, add 3 nu decay
39c5407 (nickkamp1): implement sampling for secondary decay
3d95c4e (nickkamp1): more EW decay logic, add prototype surface detector
5031575 (Nicholas Kamp): initial suite of changes to re-implement the dipole DIS from spline class in the SIREN parlance
5e2b658 (Nicholas Kamp): start MH sampling template
6f19de6 (nickkamp1): start with diff xs
72e02b4 (Nicholas Kamp): N4Bar decay fix
7c8a8f1 (Nicholas Kamp): continue implementing diff dec width, start w sampling
8283177 (Nicholas Kamp): start implementing 3 body sampling logic
88c13f0 (Nicholas Kamp): decay widths now match literature
949125f (Nicholas Kamp): overhaul the three body decay sampling logic.. some build errors remain
96f16cc (Nicholas Kamp): build fixes, total decay width fixes
b6eba95 (Nicholas Kamp): total decay width for EW decay
b895e21 (Nicholas Kamp): add lots of sampling code
c3b7765 (Nicholas Kamp): fix build warnings, add more dipole hnl splines
dc0281b (Nicholas Kamp): Total decay widths implemented, add ND280 detector
ea13b30 (nickkamp1): Rename some classes, flush out the two body decay class
f6f8b46 (Nicholas Kamp): get the decay sampling working
feddaca (Nicholas Kamp): three body decays now run!!

Notes

This branch is the squashed cumulative diff of the listed source commits. Per-commit history is browsable on dev/HNL_DIS_clean.
Large .fits data files have been removed from the branch and are packaged separately.

Squashed diff for paths: - projects/interactions/private/HNLDecay.cxx - projects/interactions/public/SIREN/interactions/HNLDecay.h - projects/interactions/private/pybindings/HNLDecay.h - projects/interactions/private/ElectroweakDecay.cxx - projects/interactions/public/SIREN/interactions/ElectroweakDecay.h - projects/interactions/private/HNLDipoleDecay.cxx - projects/interactions/public/SIREN/interactions/HNLDipoleDecay.h - projects/interactions/private/pybindings/HNLDipoleDecay.h - projects/interactions/private/pybindings/interactions.cxx Source commits on dev/HNL_DIS_clean that contributed: 197bf14 (Nicholas Kamp): fixing build errors from the merge 1bbd678 (Nicholas Kamp): fix memory leak 2219939 (nickkamp1): first step at EW decay 285285f (Nicholas Kamp): kinematics fix to 2body decays 3625e81 (Nicholas Kamp): fix build errors, add 3 nu decay 39c5407 (nickkamp1): implement sampling for secondary decay 3d95c4e (nickkamp1): more EW decay logic, add prototype surface detector 5031575 (Nicholas Kamp): initial suite of changes to re-implement the dipole DIS from spline class in the SIREN parlance 5e2b658 (Nicholas Kamp): start MH sampling template 6f19de6 (nickkamp1): start with diff xs 72e02b4 (Nicholas Kamp): N4Bar decay fix 7c8a8f1 (Nicholas Kamp): continue implementing diff dec width, start w sampling 8283177 (Nicholas Kamp): start implementing 3 body sampling logic 88c13f0 (Nicholas Kamp): decay widths now match literature 949125f (Nicholas Kamp): overhaul the three body decay sampling logic.. some build errors remain 96f16cc (Nicholas Kamp): build fixes, total decay width fixes b6eba95 (Nicholas Kamp): total decay width for EW decay b895e21 (Nicholas Kamp): add lots of sampling code c3b7765 (Nicholas Kamp): fix build warnings, add more dipole hnl splines dc0281b (Nicholas Kamp): Total decay widths implemented, add ND280 detector ea13b30 (nickkamp1): Rename some classes, flush out the two body decay class f6f8b46 (Nicholas Kamp): get the decay sampling working feddaca (Nicholas Kamp): three body decays now run!! Co-authored-by: nickkamp1 <nwkamp@mit.edu>

Copilot

Pull request overview

Note

Copilot was unable to run its full agentic suite in this review.

Adds hadronic and electroweak decay support for HNL-related interactions, including new decay classes and Python bindings, as part of the stacked refactor of dev/HNL_DIS.

Changes:

Introduces ElectroweakDecay and HNLDipoleDecay classes (headers + implementations) with decay-width calculations and final-state sampling.
Extends HNLDecay to include hadronic/meson channels, 3-body decays, and QCD corrections.
Updates pybind11 module registration to expose the new/updated decay types to Python.

Reviewed changes

Copilot reviewed 9 out of 9 changed files in this pull request and generated 10 comments.

Show a summary per file

File	Description
projects/interactions/public/SIREN/interactions/HNLDipoleDecay.h	Public API + cereal serialization for the new dipole decay.
projects/interactions/private/HNLDipoleDecay.cxx	Implements dipole decay widths, signatures, and sampling.
projects/interactions/public/SIREN/interactions/HNLDecay.h	Extends HNLDecay API (hadronic/EW helpers, updated sampling signature).
projects/interactions/private/HNLDecay.cxx	Implements hadronic/EW/3-body decay widths and sampling logic.
projects/interactions/public/SIREN/interactions/ElectroweakDecay.h	Public API + cereal serialization for electroweak boson decays.
projects/interactions/private/ElectroweakDecay.cxx	Implements W/Z decay widths, signatures, and sampling.
projects/interactions/private/pybindings/interactions.cxx	Registers new interactions/decays in the Python module.
projects/interactions/private/pybindings/HNLDipoleDecay.h	Adds Python bindings for `HNLDipoleDecay`.
projects/interactions/private/pybindings/HNLDecay.h	Adds Python bindings for `HNLDecay`.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Copilot · 2026-04-28T04:36:43Z

+            double _dipole_coupling;
+            ChiralNature _nature;
+
+            archive(::cereal::make_nvp("PrimaryTypes", _primary_types));
+            archive(::cereal::make_nvp("HNLMass", _hnl_mass));
+            archive(::cereal::make_nvp("DipoleCoupling", _dipole_coupling));
+            archive(::cereal::make_nvp("ChiralNature", _nature));
+            construct(_hnl_mass, _dipole_coupling, _nature, _primary_types);


The cereal schema is inconsistent: save() archives DipoleCoupling as std::vector<double> and ChiralNature as an int, but load_and_construct() reads DipoleCoupling into a double and ChiralNature into an enum. This will fail to deserialize. Align types by loading a std::vector<double> for DipoleCoupling (or change save() to store a single double if that’s the intended schema) and load ChiralNature via an int followed by a static_cast<ChiralNature>.

Suggested change

double _dipole_coupling;

ChiralNature _nature;

archive(::cereal::make_nvp("PrimaryTypes", _primary_types));

archive(::cereal::make_nvp("HNLMass", _hnl_mass));

archive(::cereal::make_nvp("DipoleCoupling", _dipole_coupling));

archive(::cereal::make_nvp("ChiralNature", _nature));

construct(_hnl_mass, _dipole_coupling, _nature, _primary_types);

std::vector<double> _dipole_coupling;

int _nature_int;

archive(::cereal::make_nvp("PrimaryTypes", _primary_types));

archive(::cereal::make_nvp("HNLMass", _hnl_mass));

archive(::cereal::make_nvp("DipoleCoupling", _dipole_coupling));

archive(::cereal::make_nvp("ChiralNature", _nature_int));

construct(_hnl_mass, _dipole_coupling, static_cast<ChiralNature>(_nature_int), _primary_types);

Copilot · 2026-04-28T04:36:43Z

+    void load_and_construct(Archive & archive, cereal::construct<ElectroweakDecay> & construct, std::uint32_t version) {
+        if(version == 0) {
+            std::set<siren::dataclasses::ParticleType> _primary_types;
+            construct(_primary_types);
+            archive(::cereal::make_nvp("Decay", cereal::virtual_base_class<Decay>(construct.ptr())));
+        } else {
+            throw std::runtime_error("ElectroweakDecay only supports version <= 0!");
+        }
+    }


Deserialization drops PrimaryTypes: _primary_types is never loaded from the archive, and construct(_primary_types) is called with an empty set. This makes serialized objects lose configuration. Load PrimaryTypes before constructing (and keep the save()/load_and_construct() fields symmetric).

Copilot · 2026-04-28T04:36:43Z

+        int iup = 0;
+        for(auto ubar : UpAntiQuarks) {
+          for(auto d : DownQuarks) {
+            if (record.signature.secondary_types[0]==ubar && record.signature.secondary_types[1]==d) {
+              return V_CKM.at(std::make_pair(UpQuarks[iup],d)) * GammaW;
+            }
+          }
+          ++iup;
+        }


The W→qq' partial width factor is incorrect: the amplitude should scale with |V\_CKM|, but the width scales with |V\_CKM|^2 and includes a color factor of 3 relative to leptonic channels. Returning V_CKM * GammaW underestimates/overestimates and can even go negative if signs appear. Use 3 * std::pow(V_CKM, 2) * GammaW (or equivalent) and ensure consistency with how GammaW is defined.

Copilot · 2026-04-28T04:36:44Z

+  if (proxyRecord.signature.secondary_types[0] == dataclasses::ParticleType::EMinus ||
+      proxyRecord.signature.secondary_types[0] == dataclasses::ParticleType::MuMinus ||
+      proxyRecord.signature.secondary_types[0] == dataclasses::ParticleType::TauMinus ||
+      proxyRecord.signature.secondary_types[0] == dataclasses::ParticleType::EPlus ||
+      proxyRecord.signature.secondary_types[0] == dataclasses::ParticleType::MuPlus ||
+      proxyRecord.signature.secondary_types[0] == dataclasses::ParticleType::TauPlus) {
+    for(auto meson : PlusChargedMesons) {
+      proxyRecord.signature.secondary_types[1] = meson;
+      GammaMeson += TotalDecayWidthForFinalState(proxyRecord);
+    }
+  }
+  else if (proxyRecord.signature.secondary_types[0] == dataclasses::ParticleType::EPlus ||
+           proxyRecord.signature.secondary_types[0] == dataclasses::ParticleType::MuPlus ||
+           proxyRecord.signature.secondary_types[0] == dataclasses::ParticleType::TauPlus) {


The charge-selection logic is broken: the first if includes both negative and positive leptons, so the else if branch for positive leptons is unreachable. As a result, this always sums PlusChargedMesons even when the lepton is positively charged. Split the first condition to only match negative leptons (EMinus/MuMinus/TauMinus), and keep the second branch for positive leptons to sum MinusChargedMesons.

Copilot · 2026-04-28T04:36:44Z

+  rk::P4 pHNL(geom3::Vector3(record.primary_momentum[1], record.primary_momentum[2], record.primary_momentum[3]), record.primary_mass);
+  rk::Boost boost_to_lab = pHNL.labBoost();
+  rk::P4 k3_HNLRest(p3_HNLRest*geom3::Vector3(SinTheta3_HNLRest*cos(Phi3_HNLRest),SinTheta3_HNLRest*sin(Phi3_HNLRest),CosTheta3_HNLRest),m_alpha);
+  rk::P4 k4_HNLRest(p4_HNLRest*geom3::Vector3(SinTheta4_HNLRest*cos(Phi4_HNLRest),SinTheta4_HNLRest*sin(Phi4_HNLRest),CosTheta4_HNLRest),m_alpha);


The k4 four-vector is constructed with m_alpha instead of m_beta, which will give the wrong invariant mass/energy for the second charged lepton in the 3-body final state. Use m_beta for k4_HNLRest.

Suggested change

rk::P4 k4_HNLRest(p4_HNLRest*geom3::Vector3(SinTheta4_HNLRest*cos(Phi4_HNLRest),SinTheta4_HNLRest*sin(Phi4_HNLRest),CosTheta4_HNLRest),m_alpha);

rk::P4 k4_HNLRest(p4_HNLRest*geom3::Vector3(SinTheta4_HNLRest*cos(Phi4_HNLRest),SinTheta4_HNLRest*sin(Phi4_HNLRest),CosTheta4_HNLRest),m_beta);

Copilot · 2026-04-28T04:36:44Z

 friend cereal::access;
 protected:
-HNLDecay() {};
+HNLDecay() {crundec = new CRunDec();};


Manual ownership of CRunDec* without a visible destructor/copy-control risks leaks and double-frees (e.g., copies of HNLDecay will copy the raw pointer). Prefer storing CRunDec by value or using std::unique_ptr<CRunDec> and implementing/defining appropriate copy/move behavior (or deleting copy operations).

Copilot · 2026-04-28T04:36:44Z

    ChiralNature nature;
-    const std::set<siren::dataclasses::ParticleType> primary_types = {siren::dataclasses::ParticleType::NuF4, siren::dataclasses::ParticleType::NuF4Bar};
+    const std::set<siren::dataclasses::ParticleType> primary_types = {siren::dataclasses::ParticleType::N4, siren::dataclasses::ParticleType::N4Bar};
+    CRunDec * crundec;


Manual ownership of CRunDec* without a visible destructor/copy-control risks leaks and double-frees (e.g., copies of HNLDecay will copy the raw pointer). Prefer storing CRunDec by value or using std::unique_ptr<CRunDec> and implementing/defining appropriate copy/move behavior (or deleting copy operations).

Copilot · 2026-04-28T04:36:45Z

+    HNLDecay(double hnl_mass, std::vector<double> mixing, ChiralNature nature) : hnl_mass(hnl_mass), mixing(mixing), nature(nature) {crundec = new CRunDec();};
+    HNLDecay(double hnl_mass, std::vector<double> mixing, ChiralNature nature, std::set<siren::dataclasses::ParticleType> const & primary_types) : hnl_mass(hnl_mass), mixing(mixing), nature(nature), primary_types(primary_types) {crundec = new CRunDec();};
+    HNLDecay(double hnl_mass, double mixing, ChiralNature nature) : hnl_mass(hnl_mass), mixing(std::vector<double>{0,0,mixing}), nature(nature) {crundec = new CRunDec();};
+    HNLDecay(double hnl_mass, double mixing, ChiralNature nature, std::set<siren::dataclasses::ParticleType> const & primary_types) : hnl_mass(hnl_mass), mixing(std::vector<double>{0,0,mixing}), nature(nature), primary_types(primary_types) {crundec = new CRunDec();};


Manual ownership of CRunDec* without a visible destructor/copy-control risks leaks and double-frees (e.g., copies of HNLDecay will copy the raw pointer). Prefer storing CRunDec by value or using std::unique_ptr<CRunDec> and implementing/defining appropriate copy/move behavior (or deleting copy operations).

Copilot · 2026-04-28T04:36:45Z

+    else {
+      std::cout << "HNL decay signature not recongized! Exiting\n";
+      exit(0);
+    }


Calling exit(0) inside library code makes failures non-recoverable for callers (including Python bindings) and bypasses normal error handling/cleanup. Replace exit(0) with a thrown exception (e.g., std::runtime_error) or the project’s error mechanism so the caller can handle invalid signatures.

Copilot · 2026-04-28T04:36:45Z

+    geom3::UnitVector3 pHNL_dir = pHNL_mom.direction();
+    geom3::Rotation3 x_to_pHNL_rot = geom3::rotationBetween(x_dir, pHNL_dir);
+
+    double phi = random->Uniform(0, 2.0 * M_PI);


M_PI is not guaranteed to be defined by <cmath> across all platforms/toolchains unless specific macros are set. Since this project already has siren::utilities::Constants::pi, prefer using that constant (or std::numbers::pi if available) for portability.

Suggested change

double phi = random->Uniform(0, 2.0 * M_PI);

double phi = random->Uniform(0, 2.0 * siren::utilities::Constants::pi);

austinschneider · 2026-04-28T15:32:17Z

Closing to reopen from the commit author's account so they can be added as a reviewer. Branch unchanged; will be re-linked from the new PR shortly.

Copilot AI review requested due to automatic review settings April 28, 2026 04:29

austinschneider force-pushed the feat/hnl-decays-2body branch from 9730149 to 6308e14 Compare April 28, 2026 04:34

austinschneider force-pushed the feat/hnl-decays-hadronic-ew branch from ff12237 to 26c739b Compare April 28, 2026 04:34

Copilot AI reviewed Apr 28, 2026

View reviewed changes

austinschneider force-pushed the feat/hnl-decays-2body branch from 6308e14 to cb29e52 Compare April 28, 2026 05:09

austinschneider force-pushed the feat/hnl-decays-hadronic-ew branch from 26c739b to 4aa1586 Compare April 28, 2026 05:09

austinschneider force-pushed the feat/hnl-decays-2body branch from cb29e52 to dffa635 Compare April 28, 2026 06:04

austinschneider force-pushed the feat/hnl-decays-hadronic-ew branch from 4aa1586 to 770d5fc Compare April 28, 2026 06:04

austinschneider changed the base branch from feat/hnl-decays-2body to feat/hnl-dis-from-spline April 28, 2026 13:08

austinschneider mentioned this pull request Apr 28, 2026

HNL 2-body decays #109

Closed

austinschneider force-pushed the feat/hnl-dis-from-spline branch from dffa635 to ed52cca Compare April 28, 2026 13:26

austinschneider force-pushed the feat/hnl-decays-hadronic-ew branch from 770d5fc to 88fa558 Compare April 28, 2026 13:26

austinschneider force-pushed the feat/hnl-dis-from-spline branch from ed52cca to 539374b Compare April 28, 2026 14:36

austinschneider force-pushed the feat/hnl-decays-hadronic-ew branch from 88fa558 to 73a93fd Compare April 28, 2026 14:36

austinschneider closed this Apr 28, 2026

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

HNL hadronic and electroweak decays#110

HNL hadronic and electroweak decays#110
austinschneider wants to merge 1 commit into
feat/hnl-dis-from-splinefrom
feat/hnl-decays-hadronic-ew

austinschneider commented Apr 28, 2026

Uh oh!

Copilot AI left a comment

Uh oh!

Copilot AI Apr 28, 2026

Uh oh!

Copilot AI Apr 28, 2026

Uh oh!

Copilot AI Apr 28, 2026

Uh oh!

Copilot AI Apr 28, 2026

Uh oh!

Copilot AI Apr 28, 2026

Uh oh!

Copilot AI Apr 28, 2026

Uh oh!

Copilot AI Apr 28, 2026

Uh oh!

Copilot AI Apr 28, 2026

Uh oh!

Copilot AI Apr 28, 2026

Uh oh!

Copilot AI Apr 28, 2026

Uh oh!

austinschneider commented Apr 28, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

	rk::P4 k4_HNLRest(p4_HNLRestgeom3::Vector3(SinTheta4_HNLRestcos(Phi4_HNLRest),SinTheta4_HNLRest*sin(Phi4_HNLRest),CosTheta4_HNLRest),m_alpha);
	rk::P4 k4_HNLRest(p4_HNLRestgeom3::Vector3(SinTheta4_HNLRestcos(Phi4_HNLRest),SinTheta4_HNLRest*sin(Phi4_HNLRest),CosTheta4_HNLRest),m_beta);

	double phi = random->Uniform(0, 2.0 * M_PI);
	double phi = random->Uniform(0, 2.0 * siren::utilities::Constants::pi);

Uh oh!

Conversation

austinschneider commented Apr 28, 2026

Stack: PR 12 of 14

Squashed contents

Notes

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

austinschneider commented Apr 28, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants