Neutrino Solution (NuSol) Module

The NuSol sub-package provides analytical single and double neutrino reconstruction using ellipse intersection methods. It exposes both a high-level C++ nusol class and a batch-capable CUDA/LibTorch interface (nusol_ namespace — see Neutrino Reconstruction (pyc)).

Entry Point

nusol_enum selects between the two reconstruction back-ends:

enum class nusol_enum

Selects the neutrino reconstruction algorithm.

Values:

enumerator ellipse: Use the analytic ellipse intersection method (arxiv:1305.1878).

enumerator conuix: Use the generalised n-neutrino Conuix method.

enumerator undefined: No algorithm selected.

struct nusol_t

Run parameters for the neutrino reconstruction.

Public Members

double met = 0: Missing transverse energy magnitude in MeV.

double phi = 0: Azimuthal angle of the MET vector.

double mt = 172.68 * 1000: Top-quark mass constraint in MeV (default 172.68 GeV).

double mw = 80.385 * 1000: W-boson mass constraint in MeV (default 80.385 GeV).

double violation = 0.00001: Acceptable mass-constraint violation (default 1e-5).

double limit = 0.1: Maximum allowed ΔR between the MET vector and the reconstructed neutrino (default 0.1).

int iterations = 10: Maximum number of solver iterations (default 10).

nusol_enum mode = nusol_enum::undefined : Algorithm selector.

std::vector<particle_template*> *targets = nullptr: Pointer to the list of input particles (b-quarks and leptons).

std::vector<std::pair<particle_template*, particle_template*>> *phys_pairs = nullptr: Pointer to a list of pre-specified (b, lepton) pairs. If nullptr, all combinations are tried.

class nusol : private notification, private tools

Top-level neutrino reconstruction orchestrator.

Inherits notification and tools. Dispatches to ellipse or conuix based on nusol_t::mode.

Public Functions

nusol(nusol_t *parameters)

Construct the solver with the given parameters.

Parameters:: parameters – Pointer to nusol_t.

std::vector<particle_template*> solve()

Run the reconstruction and return reconstructed neutrino candidates.

Returns:: Vector of dynamically allocated particle_template objects (caller takes ownership).

~nusol(): Destructor.

Ellipse Utilities

mtx is a minimal dense matrix class used throughout the analytical ellipse-intersection solver. All matrix operations are implemented without external BLAS dependencies.

class mtx

Dense real matrix with element-wise validity tracking.

Public Functions

mtx(): Default constructor (zero-sized).

mtx(int idx, int idy)

Construct an idx × idy matrix initialised to zero.

Parameters:

idx – Number of rows.
idy – Number of columns.

mtx(mtx *in)

Copy-construct from a pointer.

Parameters:: in – Source matrix.

mtx(mtx &in)

Copy-construct from a reference.

Parameters:: in – Source matrix.

~mtx(): Destructor.

double trace()

Compute the matrix trace.

Returns:: Sum of diagonal elements.

double det()

Compute the matrix determinant.

Returns:: Determinant.

mtx copy()

Return a deep copy.

Returns:: Copy of *this.

void copy(const mtx *ipt, int idx, int idy = -1)

Copy row idx (and optionally column idy) from ipt.

Parameters:

ipt – Source.
idx – Row index.
idy – Column index (-1 = all columns).

void copy(const mtx *ipt, int idx, int jdy, int idy)

bool assign(int idx, int idy, double val, bool valid = true)

Set element (idx,idy) to val with validity flag valid.

Returns:: true on success.

bool unique(int id1, int id2, double v1, double v2)

Set a pair of elements only if they are not already equal to v1 and v2 respectively.

Returns:: true if the assignment was performed.

mtx *cat(const mtx *v2)

Concatenate v2 as new columns to *this.

Returns:: Pointer to a new matrix (caller must delete).

mtx *slice(int idx)

Extract column idx as a column vector.

Parameters:: idx – Column index.
Returns:: Pointer to a new column-vector matrix.

mtx *eigenvalues()

Compute eigenvalues of a 3×3 matrix.

Returns:: Pointer to a 3×1 matrix of eigenvalues.

mtx *eigenvector()

Compute eigenvectors of a 3×3 matrix.

Returns:: Pointer to a 3×3 matrix whose columns are eigenvectors.

mtx T()

Return the transpose.

Returns:: Transposed matrix.

mtx inv()

Return the matrix inverse.

Returns:: Inverse.

mtx cof()

Return the cofactor matrix.

Returns:: Cofactor matrix.

mtx diag()

Return a diagonal matrix constructed from the diagonal of *this.

Returns:: Diagonal matrix.

mtx dot(const mtx &other)

Matrix-matrix product.

Parameters:: other – Right-hand factor.
Returns:: Product matrix.

mtx dot(const mtx *other)

Matrix-matrix product (pointer overload).

Parameters:: other – Pointer to the right-hand factor.
Returns:: Product matrix.

mtx cross(mtx *r1)

Cross product of two 3-vectors.

Parameters:: r1 – Second vector.
Returns:: Cross-product vector.

mtx cross(mtx *r1, mtx *r2)

bool valid(int idx, int idy)

Test whether element (idx,idy) has a valid flag.

Returns:: true if valid.

void print(int prec = 6, int width = 12)

Print the matrix to stdout with prec decimal digits.

Parameters:

prec – Decimal precision (default 6).
width – Column width (default 12).

mtx &operator=(const mtx &o)

double m_00()

double m_01()

double m_02()

double m_10()

double m_11()

double m_12()

double m_20()

double m_21()

double m_22()

Public Members

double **_m = nullptr: Raw 2D array of matrix elements.

bool **_b = nullptr: Validity mask (same dimensions as _m).

int dim_i = 0: Number of rows.

int dim_j = 0: Number of columns.

double tol = 1e-9: Numerical tolerance for near-zero tests (default 1e-9).

Constrained Neutrino Solutions

conuix orchestrates multi-neutrino reconstruction by creating one conuic instance per lepton–b-jet pair and merging the solutions.

class conuix

Orchestrates multi-pair Conuix neutrino reconstruction.

Public Functions

conuix(nusol_t *params)

Construct and initialise all conuic sub-solvers.

Parameters:: params – Pointer to the nusol_t run parameters.

~conuix(): Destructor. Frees all owned conuic objects.

std::vector<particle_template*> nunu_make()

Run all sub-solvers and build the reconstructed neutrino candidates.

Returns:: Vector of reconstructed particle_template objects.

Public Members

std::vector<conuic*> *cnx = nullptr: List of conuic solvers (one per b–lepton pair).

nusol_t *params = nullptr: Pointer to the run parameters.

std::string prefix = "": Optional log prefix string.

conuic implements the single-neutrino constrained solver. It finds the neutrino momentum that simultaneously satisfies the W-mass and top-mass constraints using the Möbius-parametrisation approach.

class conuic

Single b–lepton neutrino solver via Möbius / pencil-of-conics.

Public Functions

conuic(particle_template *jet, particle_template *lep)

Construct the solver for the (b-jet, lepton) pair.

Parameters:

jet – Pointer to the b-jet.
lep – Pointer to the lepton.

~conuic(): Destructor.

long double Z2(long double Sx, long double Sy)

Evaluate Z² from the shift parameters.

Parameters:

Sx – x-shift.
Sy – y-shift.

Returns:

Z².

long double Sx(long double Tau, long double Z)

Evaluate S_x from (τ, Z).

Parameters:

Tau – τ parameter.
Z – Z parameter.

Returns:

S_x.

long double Sy(long double Tau, long double Z)

Evaluate S_y from (τ, Z).

Returns:: S_y.

long double x1(long double t, long double Z): Return the x₁ coordinate of the neutrino.

long double y1(long double t, long double Z): Return the y₁ coordinate of the neutrino.

bool get_TauZ(long double sx, long double sy, long double *z, long double *t)

Compute (τ, Z) from the given (S_x, S_y) values.

Parameters:

sx – S_x value.
sy – S_y value.
z – Output Z.
t – Output τ.

Returns:

true on success.

long double P(long double l, long double t, long double Z)

Evaluate the characteristic polynomial P(λ, τ, Z).

Parameters:

l – Lambda.
t – τ.
Z – Z.

Returns:

P value.

long double dPdt(long double l, long double t, long double Z)

Evaluate ∂P/∂τ.

Returns:: Derivative value.

long double dPdtL0(long double t, long double Z)

long double dPl0(long double t)

void dPdZ0(long double z, long double t, std::complex<long double> *l1, std::complex<long double> *l2)

matrix_t Hmatrix(long double t, long double Z)

matrix_t Nmatrix(long double t, long double Z)

bool mass_line(long double mW, long double mT, std::complex<long double> *z_ = nullptr, long double *t_ = nullptr)

Find the point on the mass line (mW, mT) closest to the neutrino ellipse.

Parameters:

mW – W-boson mass constraint.
mT – Top-quark mass constraint.
z_ – Output Z value.
t_ – Output τ value.

Returns:

true if a solution was found.

void debug(): Print diagnostic information to stdout.

Public Members

long double tstar = 0: Optimal τ value found by the solver.

long double error = 0: Residual penalty at tstar.

long double theta = 0: Angle associated with the optimal solution.

bool converged = false: true if the solver converged.

matrix_t vstar: Neutrino four-momentum vector at the optimal tstar.

atomics_t *cache = nullptr: Pre-computed kinematic constants for this pair.

Conuix Internal Structs (`Conuix` namespace)

struct atomics_t : public Conuix::debug

Master cache struct that pre-computes all kinematic constants for a single (b-jet, lepton) pair needed by the Conuix solver.

Inherits Conuix::debug.

Public Functions

atomics_t(particle_template *jet, particle_template *lep, double m_nu = 0)

Construct and pre-compute all constants.

Parameters:

jet – Pointer to the b-jet.
lep – Pointer to the lepton.
m_nu – Neutrino mass hypothesis (default 0).

~atomics_t() override: Destructor.

long double x1(long double Z, long double t): Evaluate the x₁ neutrino coordinate.

long double y1(long double Z, long double t): Evaluate the y₁ neutrino coordinate.

bool GetTauZ(long double sx, long double sy, long double *z, long double *t)

Invert (S_x, S_y) → (τ, Z).

Returns:: true on success.

void masses(long double Z, long double t, std::complex<long double> *mw, std::complex<long double> *mt)

Compute the W and top masses corresponding to (Z, τ).

Parameters:

Z – Z parameter.
t – τ parameter.
mw – Output W mass.
mt – Output top mass.

void eigenvector(long double tau, matrix_t *vstar, long double *theta_s)

Compute the eigenvector at tau and its associated angle.

Parameters:

tau – τ value.
vstar – Output eigenvector.
theta_s – Output angle.

void debug_mode(particle_template *jet, particle_template *lep)

Public Members

Conuix::kinematic_t lp

Conuix::kinematic_t jt

Conuix::kinematic_t nu

Conuix::base_t base

Conuix::rotation_t rotation

Conuix::thetapsi_t psi_theta

Conuix::pencil_t pencil

Conuix::Sx_t Sx

Conuix::Sy_t Sy

Conuix::H_matrix_t H_Matrix

Conuix::Mobius_t Mobius

struct kinematic_t : public Conuix::debug

Compact kinematics for one particle (β, mass, energy).

Public Functions

virtual void print(int p = 16) override

Public Members

long double beta = 0

long double mass = 0

long double energy = 0

struct rotation_t : public Conuix::debug

Rotation from the decay frame to the lab frame (φ, θ, R^T).

Public Functions

virtual void print(int p = 16) override

Public Members

long double phi = 0

long double theta = 0

matrix_t vec

matrix_t R_T

struct base_t : public Conuix::debug

Base kinematic constants used to construct the pencil of conics.

Stores the (ω, Ω, β, cos/sin, …) constants shared by all subsequent pencil evaluations.

Public Functions

virtual void print(int p = 16) override

Public Members

long double rbl = 0

long double cos = 0

long double sin = 0

long double w = 0

long double o = 0

long double w2 = 0

long double o2 = 0

long double beta = 0

long double mass = 0

long double E = 0

long double tpsi = 0

long double cpsi = 0

long double spsi = 0

struct pencil_t : public Conuix::debug

Pencil-of-conics polynomial coefficients.

Public Functions

long double Z2(long double Sx, long double Sy)

virtual void print(int p) override

Public Members

long double a = 0

long double b = 0

long double c = 0

long double d = 0

long double e = 0

struct Sx_t : public Conuix::debug

Polynomial coefficients for the S_x shift.

Public Functions

long double Sx(long double tau, long double Z)

virtual void print(int p) override

Public Members

long double a = 0

long double b = 0

long double c = 0

struct Sy_t : public Conuix::debug

Polynomial coefficients for the S_y shift.

Public Functions

long double Sy(long double tau, long double Z)

virtual void print(int p) override

Public Members

long double a = 0

long double b = 0

long double c = 0

struct H_matrix_t : public Conuix::debug

H̄ and H matrix components (tau-independent, sin, cos parts).

Public Functions

matrix_t H_Matrix(long double tau, long double Z)

matrix_t H_Tilde(long double tau, long double Z)

virtual void print(int p) override

Public Members

matrix_t HBX

matrix_t HBS

matrix_t HBC

matrix_t HTX

matrix_t HTS

matrix_t HTC

struct Mobius_t : public Conuix::debug

Möbius transform coefficients and the τ* minimiser.

The Möbius transform maps the conic pencil parameter λ to the neutrino parameter τ. Mobius_t stores the four Möbius coefficients (a, b, c, d), the eigenvalues, and the result of the τ* minimisation.

Public Functions

void init(base_t *base)

long double alpha_p(long double u)

long double alpha_m(long double u)

long double P(long double z, long double l, long double t)

long double dPdt(long double z, long double l, long double t)

long double dPdtL0(long double z, long double t)

long double dPl0(long double x, bool use_tanh = false)

Public Members

long double a

long double b

long double c

long double d

long double pole1 = 0

long double pole2 = 0

long double a_

long double b_

std::complex<long double> eig_v1

std::complex<long double> eig_v2

std::complex<long double> eig_vx[2]

std::complex<long double> eig_vy[2]

std::complex<long double> kfactor

std::complex<long double> f1_pts

std::complex<long double> f2_pts

std::complex<long double> mid

long double error = 0

long double tstar = 0

bool converged = false

ODE Runge–Kutta Multi-Neutrino Solver

odeRK implements a fourth-order Runge–Kutta ODE integrator used by the multi-solution multisol back-end to evolve ellipse states until convergence.

class odeRK : public tools

RK4 integrator for multi-neutrino ellipse alignment.

Inherits tools.

Public Functions

odeRK(std::vector<multisol*> *sols, vec3 met, int iter, double step)

Construct the integrator.

Parameters:

sols – Pointer to the list of multisol surfaces to integrate.
met – MET three-vector.
iter – Maximum number of RK4 steps.
step – Step size Δt.

~odeRK(): Destructor.

void solve(): Run the full RK4 integration and populate _state.

void rk4(double dt)

Perform one RK4 step with step size dt.

Parameters:: dt – Step size.

void update_t(): Update all t parameters from the current state.

double solve_z_phi()

Solve for the z and φ parameters that minimise the MET residual.

Returns:: Residual at the solution.

double residual(std::vector<double> wg, std::vector<double> phx)

Evaluate the MET residual for the given angles.

Parameters:

wg – Weights / z values.
phx – φ values.

Returns:

Residual.

std::vector<ellipse_t> derivative(const std::vector<ellipse_t> &dS)

Compute the time derivative of the state dS.

Returns:: Derivative state.

double ghost_angle(int nui)

std::vector<double> plane_rk4(const std::vector<double> &t_initial)

std::vector<double> plane_align(const std::vector<ellipse_t> &current_state)

struct ellipse_t

State vector for one neutrino ellipse at a given (t, z).

Public Functions

void print()

Public Members

vec3 A: Ellipse basis vector A.

vec3 B: Ellipse basis vector B.

vec3 C: Ellipse centre vector C.

vec3 vA

vec3 vB

vec3 vC

double t: Current t parameter.

double z = 1.0: Current z parameter (default 1.0).

struct recon_t

Reconstruction result for a single di-neutrino configuration.

Public Members

bool is_valid = false: true if a valid solution was found.

double residual = -1: MET residual (lower = better; default -1 = not computed).

std::vector<double> t: Optimal t values for each neutrino.

std::vector<double> z: Optimal z values for each neutrino.

std::vector<double> phi: Optimal φ values for each neutrino.

Neutrino Solution (NuSol) Module

Entry Point

Ellipse Utilities

Constrained Neutrino Solutions

Conuix Internal Structs (Conuix namespace)

ODE Runge–Kutta Multi-Neutrino Solver

Conuix Internal Structs (`Conuix` namespace)