Consolidate nudk decay-osc calculator (supersedes #42, #47) by gavinsdavies · Pull Request #48 · cafana/OscLib

gavinsdavies · 2026-07-29T19:32:14Z

Andrea (acbarros) opened two PRs for the same nudk decay-oscillation calculator work — #42 (with the full review history from @gavinsdavies and @cullenmsullivan) and #47 (a later rebuild off main that addressed the using namespace std; / style feedback and fixed fIsNuBar handling, but lost the review thread and had zero reviews of its own).

This branch takes #42 as the base (preserving the review history) and applies #47's fixes on top as a single new commit, so reviewers aren't starting over. One typo introduced in #47 ("Neutirno" -> "Neutrino") was fixed along the way.

Pushed to a branch on cafana/OscLib rather than force-pushing onto Andrea's personal fork, to avoid rewriting their fork history.

Next: will work through @cullenmsullivan's remaining style comments from #42 here to help Andrea get this over the line. Closes #42 and #47 in favor of this.

Copy Cullen changes

…after all comments Merge branch 'main' into feature/acbarros_nudk_osc

Merge branch 'feature/acbarros_nudk_osc' of https://github.com/abarroscs/OscLib into feature/acbarros_nudk_osc

Consolidates the fixes from PR #47 onto this branch so the existing review thread stays attached to the current code. Fixes a typo introduced in #47 (Neutirno -> Neutrino). Pushed to a branch on cafana/OscLib rather than force-pushing onto Andrea's personal fork (abarroscs/OscLib) to avoid rewriting their fork history.

Strips the commented-out three-argument SolveHam overload (.cxx/.h) and its matching //int anti; member, per cullenmsullivan's review on #42. Also drops the now-stale @param anti doxygen line that documented the removed overload.

gavinsdavies · 2026-07-29T19:40:32Z

Went through @cullenmsullivan's review comments from #42 against the current code:

assert(false) → abort() — already done
Missing email for Mario Acero — already filled in
Inconsistent cout/std::cout, drop using namespace std; — fixed by the consolidation commit
Two "crash on user error" threads — already marked resolved, matches current code
Stale fBuiltHms = false; comment — already removed in an earlier This are nudk codes #42 commit
Commented-out SolveHam(double E, double Ne, int anti) — removed just now, along with the matching declaration/member in the header and the now-stale @param anti doxygen line

All of @cullenmsullivan's outstanding points on the old thread are addressed. Remaining open items are the cosmetic ones from my own #42 review (naming, stray includes, an incorrect class description) — will follow up on those separately.

The class brief was updated to say 'decay oscillation calculators' but the explanatory paragraph below it still described a 'sterile oscillation calculator', copy-pasted from IOscCalcSterile.h.

gavinsdavies · 2026-07-29T19:46:13Z

Went through my own #42 review comments as well:

Stray absolute-path includes (/exp/nova/app/users/acbarros/...) in IOscCalcDecay.cxx/.h and OscCalcDecayEigen.h — already removed
Dead commented-out SetAlpha3/SetAlpha2 overloads in IOscCalcDecay.cxx — already removed
calc_sterile → calc_decay naming in DowncastToDecay — already fixed
Dead commented GetAlpha3()/GetAlpha2() and SetIsNuBar(int)/GetIsNuBar() declarations in OscCalcDecayEigen.h — already removed
fGotES/kr2GNe question — fGotES now has a doc comment, kr2GNe was dropped (unused)
Fixed just now: the IOscCalcDecay class doc comment's \brief line said "decay oscillation calculators" but the paragraph underneath still described a "sterile oscillation calculator" (leftover from copy-pasting IOscCalcSterile.h's doc block) — reworded to match.

All review comments from #42 (mine and @cullenmsullivan's) are now addressed on this branch. Should be ready for a fresh look.

They duplicate storage and getters already provided by IOscCalcAdjustable.

chenel

I had Claude take a look at this PR also before signing off. It found a few more issues, below.
There are 3 things below (fully removing fGotES, fully removing the unused fEval/eigen-solve state, and the dead 2-flavor branch) span several non-adjacent sites, so the suggestion covers the primary edit and the comment describes the rest.

Items 1 and 2 are the most important ones before merging, since they concern uninitialized behavior, but since it's almost all just confirming that Claude's suggestions are ok and clicking "apply suggestion", hopefully it's all straightforward.

Co-authored-by: Jeremy Wolcott <chenel@users.noreply.github.com>

@chenel

complexsolver() filled fEval with the Hamiltonian's eigenvalues, but fEval was never read anywhere (propagation uses fHam.exp() in PropagatePath instead), so every uncached energy point paid for a ComplexEigenSolver call whose result was discarded. The paired fGotES cache flag was never checked either. Removes both members, all their writers (constructors, SetIsNuBar, BuildHms), the complexsolver method itself, and the now-dead call site in SolveHam. Per @chenel's review on #48.

cullenmsullivan

Looks mostly good to me - my primary suggestion is to avoid using == to check for equivalence of floating point variables. I've suggested a different method that's safer.

cullenmsullivan · 2026-08-01T15:48:49Z

+  }
+
+  //---------------------------------------------------------------------------    
+  const IOscCalcDecay* DowncastToDecay(const IOscCalc* calc, bool quiet)


Suggested change

const IOscCalcDecay* DowncastToDecay(const IOscCalc* calc, bool quiet)

const IOscCalcDecay* OscCalcDecayTrivial::DowncastToDecay(const IOscCalc* calc, bool quiet)

(However, it's unclear to me what the use case is, and it might be worth asking if this is even needed or whether it should belong in OscCalcDecay instead of OscCalcDecayTrivial is.)

cullenmsullivan · 2026-08-01T15:49:43Z

+    return nullptr; // If the cast failed, calc_decay should be nullptr anyway
+  }
+
+  IOscCalcDecay* DowncastToDecay(IOscCalc* calc, bool quiet)


Suggested change

IOscCalcDecay* DowncastToDecay(IOscCalc* calc, bool quiet)

IOscCalcDecay* OscCalcTrivial::DowncastToDecay(IOscCalc* calc, bool quiet)

cullenmsullivan · 2026-08-01T16:04:16Z

+      this->SetL(810);
+      this->SetRho(2.74);
+    }
+    else if(fNumNus==2){


This condition will never be met because, as far as I can tell, fNumNus is always initialized to 3, and there's no setter. Suggest either removing the condition or adding a setter method.

cullenmsullivan · 2026-08-01T16:04:55Z

+    if (i < 1 || i > fNumNus - 1 || j < 2 || j > fNumNus) {
+      std::cout << "Fatal Error Occurred" << std::endl;
+      std::cout << "Theta" << i << j << " not valid for " << fNumNus;
+      std::cout << " neutrinos. Abortig." << std::endl;


Suggested change

std::cout << " neutrinos. Abortig." << std::endl;

std::cout << " neutrinos. Aborting." << std::endl;

Pedantic: fixing a typo

cullenmsullivan · 2026-08-01T16:13:08Z

+      abort();
+    }
+    // Check if value is actually changing                                                                                                                             
+    fBuiltHms = fBuiltHms && (fTheta(i - 1,j - 1) == th);


Suggested change

fBuiltHms = fBuiltHms && (fTheta(i - 1,j - 1) == th);

fBuiltHms = fBuiltHms && (std::abs(fTheta(i - 1,j - 1) - th) <= std::numeric_limits<double>::epsilon());

Does this work? I thiiiiink this is unsafe - you're checking for equivalence between two floating point values using raw ==. See: https://stackoverflow.com/questions/4548004/how-to-correctly-and-standardly-compare-floats

cullenmsullivan · 2026-08-01T16:14:18Z

+      abort();
+    }
+    // Check if value is actually changing                                                                                                                             
+    fBuiltHms = fBuiltHms && (fDelta(i - 1,j - 1) == delta);


Suggested change

fBuiltHms = fBuiltHms && (fDelta(i - 1,j - 1) == delta);

fBuiltHms = fBuiltHms && (std::abs(fDelta(i - 1,j - 1) - delta) <= std::numeric_limits<double>::epsilon());

Another floating point equivalence check written safer.

cullenmsullivan · 2026-08-01T16:14:46Z

+      return;
+    }
+    // Check if value is actually changing                                                                                                                             
+    fBuiltHms = fBuiltHms && (fDm(i - 1) == dm);


Suggested change

fBuiltHms = fBuiltHms && (fDm(i - 1) == dm);

fBuiltHms = fBuiltHms && (std::abs(fDm(i - 1) - dm) <= std::numeric_limits<double>::epsilon());

cullenmsullivan · 2026-08-01T16:15:44Z

+  void OscCalcDecayEigen::RotateH(int i, int j, Eigen::Matrix3cd& Ham)
+  {
+    // Do nothing if angle is zero                                                                                                                                     
+    if (fTheta(i,j) == 0) return;


Suggested change

if (fTheta(i,j) == 0) return;

if (std::abs(fTheta(i,j)) <= std::numeric_limits<double>::epsilon()) return;

cullenmsullivan · 2026-08-01T16:17:13Z

+
+  void OscCalcDecayEigen::SolveHam(double E, double Ne)
+  {
+    if(Ne!=fCachedNe || E!=fCachedE || !fBuiltHms ){


Suggested change

if(Ne!=fCachedNe || E!=fCachedE || !fBuiltHms ){

if(std::abs(Ne-fCachedNe) <= std::numeric_limits<double>::epsilon() || std::abs(E-fCachedE) <= std::numeric_limits<double>::epsilon() || !fBuiltHms ){

cullenmsullivan · 2026-08-01T16:19:57Z

+#include <iostream>
+#include <cassert>
+#include <stdlib.h>
+#include <complex>


Suggested change

#include <complex>

#include <complex>

#include <cmath>

#include <limits>

These headers are used for the safer floating point equivalence checks.

Andrea Connelly and others added 12 commits March 2, 2026 08:53

This are nudk codes

0d88b82

This is Decay Calc

e1d639c

This is the update for NuDecay

13a57c0

Merge pull request #1 from cafana/main

bee6229

Copy Cullen changes

Accepting the changes made for Cullen and adding the nudecay updates …

05dc276

…after all comments Merge branch 'main' into feature/acbarros_nudk_osc

changes in lines 108 128 185 212 227

504b887

Merge branch 'feature/acbarros_nudk_osc' of https://github.com/abarroscs/OscLib into feature/acbarros_nudk_osc

changes made OscCalcDecayEigen.cxx

a2a8374

Merge branch 'feature/acbarros_nudk_osc' of https://github.com/abarroscs/OscLib into feature/acbarros_nudk_osc

Test comment

9d932e2

Update IOscCalcDecay.h

5f773e9

Update OscCalcDecayEigen.cxx

19bf681

Update OscCalcDecayEigen.h

bab4124

gavinsdavies requested a review from a team as a code owner July 29, 2026 19:32

This was referenced Jul 29, 2026

This are nudk codes #42

Closed

Changes made in OscCalcDecayEigen.cxx #47

Closed

Remove dead SolveHam(..., int anti) declarations

5661c36

Strips the commented-out three-argument SolveHam overload (.cxx/.h) and its matching //int anti; member, per cullenmsullivan's review on #42. Also drops the now-stale @param anti doxygen line that documented the removed overload.

Fix stale sterile-calculator wording in IOscCalcDecay doc comment

9d27665

The class brief was updated to say 'decay oscillation calculators' but the explanatory paragraph below it still described a 'sterile oscillation calculator', copy-pasted from IOscCalcSterile.h.

cullenmsullivan reviewed Jul 29, 2026

View reviewed changes

Comment thread OscLib/IOscCalcDecay.h Outdated

Drop redundant fRho/GetL/GetRho from IOscCalcDecay

a7198cc

They duplicate storage and getters already provided by IOscCalcAdjustable.

This comment was marked as outdated.

Sign in to view

chenel reviewed Jul 30, 2026

View reviewed changes

gavinsdavies and others added 7 commits July 31, 2026 10:10

Update OscLib/OscCalcDecayEigen.cxx

2af34ab

Co-authored-by: Jeremy Wolcott <chenel@users.noreply.github.com>

Update OscLib/OscCalcDecayEigen.cxx

1cfddaa

Co-authored-by: Jeremy Wolcott <chenel@users.noreply.github.com>

Update OscLib/IOscCalcDecay.h

5558123

Co-authored-by: Jeremy Wolcott <chenel@users.noreply.github.com>

Update OscLib/IOscCalcDecay.h

315bedb

Co-authored-by: Jeremy Wolcott <chenel@users.noreply.github.com>

Update OscLib/OscCalcDecayEigen.cxx

2c86c5f

Co-authored-by: Jeremy Wolcott <chenel@users.noreply.github.com>

Update OscLib/OscCalcDecayEigen.cxx

7560973

Co-authored-by: Jeremy Wolcott <chenel@users.noreply.github.com>

Update OscLib/OscCalcDecayEigen.h

0734ce0

Co-authored-by: Jeremy Wolcott <chenel@users.noreply.github.com>

gavinsdavies and others added 2 commits July 31, 2026 10:13

Update OscLib/OscCalcDecayEigen.h

7e1f983

Co-authored-by: Jeremy Wolcott <chenel@users.noreply.github.com>

Update OscLib/OscCalcDecayEigen.cxx

d43b223

Co-authored-by: Jeremy Wolcott <chenel@users.noreply.github.com>

gavinsdavies mentioned this pull request Jul 31, 2026

Remove dead fEval/fGotES eigenvalue-solve machinery #49

Open

Update OscLib/OscCalcDecayEigen.h

ac78ddc

Co-authored-by: Jeremy Wolcott <chenel@users.noreply.github.com>

gavinsdavies requested a review from cullenmsullivan July 31, 2026 16:25

cullenmsullivan suggested changes Aug 1, 2026

View reviewed changes

	const IOscCalcDecay* DowncastToDecay(const IOscCalc* calc, bool quiet)
	const IOscCalcDecay* OscCalcDecayTrivial::DowncastToDecay(const IOscCalc* calc, bool quiet)

	IOscCalcDecay* DowncastToDecay(IOscCalc* calc, bool quiet)
	IOscCalcDecay* OscCalcTrivial::DowncastToDecay(IOscCalc* calc, bool quiet)

	std::cout << " neutrinos. Abortig." << std::endl;
	std::cout << " neutrinos. Aborting." << std::endl;

	fBuiltHms = fBuiltHms && (fTheta(i - 1,j - 1) == th);
	fBuiltHms = fBuiltHms && (std::abs(fTheta(i - 1,j - 1) - th) <= std::numeric_limits<double>::epsilon());

	fBuiltHms = fBuiltHms && (fDelta(i - 1,j - 1) == delta);
	fBuiltHms = fBuiltHms && (std::abs(fDelta(i - 1,j - 1) - delta) <= std::numeric_limits<double>::epsilon());

	fBuiltHms = fBuiltHms && (fDm(i - 1) == dm);
	fBuiltHms = fBuiltHms && (std::abs(fDm(i - 1) - dm) <= std::numeric_limits<double>::epsilon());

	if (fTheta(i,j) == 0) return;
	if (std::abs(fTheta(i,j)) <= std::numeric_limits<double>::epsilon()) return;

	if(Ne!=fCachedNe \|\| E!=fCachedE \|\| !fBuiltHms ){
	if(std::abs(Ne-fCachedNe) <= std::numeric_limits<double>::epsilon() \|\| std::abs(E-fCachedE) <= std::numeric_limits<double>::epsilon() \|\| !fBuiltHms ){

Uh oh!

Conversation

gavinsdavies commented Jul 29, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

gavinsdavies commented Jul 29, 2026

Uh oh!

gavinsdavies commented Jul 29, 2026

Uh oh!

Uh oh!

This comment was marked as outdated.

Uh oh!

chenel left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

cullenmsullivan left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

gavinsdavies commented Jul 29, 2026 •

edited

Loading

chenel left a comment •

edited

Loading