Charm DIS injection#74
Conversation
nickkamp1
left a comment
There was a problem hiding this comment.
looks good, just need to remove some debug printouts and a few files that are used for analysis of d meson events. I will merge after testing this locally, but @austinschneider feel free to take a look and leave comments if you would like
|
Made changes! made request of new round of review, please check if the changes are made |
|
I'm starting to take a look now. I'll provide detailed comments sometime today :) |
|
It is not clear to me that we need this In the way that it has been implemented here, can't the hadronization process just be implemented as a new decay process with |
|
Regardless of if we keep the |
| double const & total_decay_length) const { | ||
| Vector3D direction = p0 - intersections.position; | ||
| if(direction.magnitude() == 0) { | ||
| if(direction.magnitude() == 0 || direction.magnitude() <= 1e-5) { |
There was a problem hiding this comment.
| if(direction.magnitude() == 0 || direction.magnitude() <= 1e-5) { | |
| // Check that direction is approximately zero | |
| if(direction.magnitude() <= 1e-5) { |
There was a problem hiding this comment.
Was there a specific issue that the distance.magnitude() <= 1e-5 checks fixed?
I'm wondering if there is another underlying issue that this is covering up, or if we really need this fix.
There was a problem hiding this comment.
Instead of adding special login in this class checks for a hadronization interaction to change the class behavior, you should implement a new secondary vertex distribution that has the appropriate behavior for hadronization.
There was a problem hiding this comment.
@MiaochenJin are we still using Secondary{Bounded,Physical}VertexDistribution for charms? I think it's probably ok as long as we add a check that any particle with hadronization interactions has no scattering or decay interactions present, since hadronization would override those
There was a problem hiding this comment.
You are right. we do not currently have any hadronization logic any more. A new commit will remove all dead code, and we shouldn't need to change the Secondary*Distribution definitions any more.
| #include "SIREN/interactions/CrossSection.h" | ||
| #include "SIREN/interactions/Decay.h" | ||
| #include "SIREN/interactions/Hadronization.h" | ||
|
|
| #include "SIREN/interactions/CrossSection.h" | ||
| #include "SIREN/interactions/Decay.h" | ||
| #include "SIREN/interactions/Hadronization.h" | ||
|
|
| #include "SIREN/interactions/CrossSection.h" | ||
| #include "SIREN/interactions/Decay.h" | ||
| #include "SIREN/interactions/Hadronization.h" | ||
|
|
| #include "SIREN/interactions/CrossSection.h" | ||
| #include "SIREN/interactions/Decay.h" | ||
| #include "SIREN/interactions/Hadronization.h" | ||
|
|
…avoid double counting
…avoid double counting
… the kinematics should be energy and momentum conserving, so the Q^2 distortion would be the only effect
The post-MH hadronic remnant p4X is built via P4 arithmetic (operator+/-) that resets the mass cache, forcing P4::m() to recompute msq = e^2 - p.lengthSquared() independently of dot(). The do-while above already accepts p4X based on dot() >= 0, but roundoff between dot() and lengthSquared() can be ~1e-15 in either direction, occasionally triggering the assert msq >= 0 inside m(). Fix: use the validated dot() result directly when setting the hadronic remnant mass, clipping at 0 for safety. The do-while remains the gatekeeper for accepting p4X. The new (xi, y) sampling reaches this edge case more frequently than the old (x, y) sampling did; observed crash at event 1569 of a deterministic 10k run with seed 1234.
Initialize pqy_lab to quiet NaN before the pqx_lab>momq_lab scaling loop and throw InjectionFailure if pqy_lab is still NaN after the loop. The loop body has multiple early-exit branches (break statements) that can leave pqy_lab unwritten if iteration_break/iteration_attempt / E1_lab relationships pin the value below threshold without finding a fix; the old code then read uninitialized memory in the subsequent Eq_lab and final-state momentum construction. Convert the silent footgun into a clean exception that the injector can catch and reject the event.
…rride Header carried an override declaration that was a leftover from an older internal version of the class where the spline was treated as inclusive and the override returned it raw. With the per-signature fragmentation fraction now applied inside TotalCrossSection (Andys 3341768), the correct behavior is the base class default: loop over the D-meson signatures and sum their per-signature TotalCrossSection. That gives spline x sum(fragfracs) = ~0.98 x spline (D0 + D+/- + Ds+/- = 0.6 + 0.23 + 0.15), with the small remainder corresponding to Lambda_c production that is not yet modeled. Dropping the override avoids a link-time unresolved symbol and lets the base default do the right thing.
Extend three helpers in QuarkDISFromSpline so Ds mesons are produced alongside D0 and D+/-: * FragmentationFraction: Ds+/Ds- = 0.15, matching the 0.60:0.23:0.15 thesis convention used in PythiaDISCrossSection and CharmHadronization. Sum across the three modeled D-types is 0.98; the residual ~0.02 corresponds to Lambda_c production not yet modeled. * DTypesForPrimary: add DsPlus to the nu set and DsMinus to the nubar set so InitializeSignatures generates symmetric c/c-bar Ds signatures. * getHadronMass: 1.96834 GeV (PDG 2022) for both DsPlus and DsMinus. Constants.h does not yet expose DsMass; literal is consistent with the Ds mass used in CharmMesonDecay::particleMass. isD() in dataclasses/Particle.cxx already recognizes Ds+/Ds- (added earlier in this branch), so the loop in TotalCrossSection that applies fragmentation fraction to D-meson signatures now picks up Ds without further changes.
After Andys 1110a2d (D0Mass / DPlusMass literal swap), Constants::D0Mass now correctly returns 1.86484 GeV (PDG D0). Three places in this branch hardcoded the previous wrong value 1.86962 (which actually corresponded to the D+ mass) for Constants::D0Mass: * QuarkDISFromSpline.h:42 - documentation comment * tests/slow_rescaling/smoke_quarkdis_100.py:13 - sampling constant * tests/slow_rescaling/smoke_quarkdis_10k.py:21 - sampling constant Update all three to 1.86484 so they match the post-fix Constants::D0Mass that the C++ side emits.
Andys upstream 3341768 removed the vestigial int quark_type parameter from all QuarkDISFromSpline ctor signatures. Our smoke tests still passed it (predating that cleanup), causing TypeError on construction. Drop the line.
* QuarkDISFromSpline::SampleFinalState: bump precision-loop max
iterations from 10 to 15. Empirically, 10 vs 15 produces the same
6/10000 rejection rate at this E_nu scale, but giving the loop a few
more chances at marginal-kinematics events is cheap insurance.
* tests/slow_rescaling/smoke_quarkdis_{100,10k}.py: on InjectionFailure
from the NaN guard, resample with fresh RNG state and retry up to
100 times before treating the event slot as unrecoverable. Mirrors
the SIREN Injector framework which retries InjectionFailure events
automatically; smoke tests bypass the Injector and need to do it
explicitly.
Track and report:
- total resample count (events that needed >= 1 retry)
- unrecoverable count (retry budget exhausted)
Result: 10k smoke OK 10000/10000 sampled with 6 resamples; 100 OK
100/100. No spurious test failures.
…o primaries
InitializeSignatures unconditionally registered {D0, D+, Ds+} for both
neutrino and antineutrino primaries. SampleFinalState mutates the
signature's meson slot to whatever Pythia actually produced — for ν̄
charm-CC that is {D̄0, D-, Ds-}, which then did not match any registered
signature. The weighter could not look the post-mutation signature back
up and event_weight came out NaN on essentially every ν̄ event.
Pick the CP-mirror set when the primary is NuEBar / NuMuBar / NuTauBar
so registered signatures stay consistent with what SampleFinalState
writes back. ν side unchanged.
Verified: with the fix, NuMuBar Phase-0 produces 100% finite weights
(20/20 in the smoke test, vs ~0/1000 before).
…metry, hardening - QuarkDISFromSpline now samples directly in (xi, y); legacy outgoing-xi remap removed. Kinematic check rebuilt around slow-rescaling identities. Precision loop hardened (NaN guard, maxiter 10->15) with InjectionFailure on non-convergence. - CharmMesonDecay rewritten to Pythia-style 3-body phase space with V-A ME accept-reject for D+/D-/D0/D0Bar; pure phase space for Ds with eta/eta-prime/phi as hadron daughter. All c-bar mirror cases added. - PythiaDISCrossSection: per-event isoscalar target (10/18 p, 8/18 n for H2O); per-event Pythia reinit so beam energy tracks E_nu; const_cast PID override writes Pythia's actually-produced meson into the signature so c-bar gets D0Bar/D-/Ds-. - Ds+/Ds- added to isD(), DTypesForPrimary, IsCharmedHadron, getIndices, FragmentationFraction (0.15), Particle::isD, CharmMesonDecay primary set, DMesonELoss primary set, CharmMesonDecay decay-mode coverage. - ParticleTypes.def: DsMinusBar -> DsMinus typo fix. - Weighter.tcc: switch totals to TotalCrossSectionAllFinalStates (no-op refactor under our Bug Harvard-Neutrino#6 fragfrac semantics; both paths return 0.98sig per the per-signature sum 0.60 + 0.23 + 0.15). - New cmake/photospline_patches/0002-move-cholmod-include-outside-extern-c for suite-sparse 7.x compatibility. - New tests/slow_rescaling/smoke_quarkdis_{100,10k}.py. Our 2026-03-24 Bug Harvard-Neutrino#6 fragfrac multiplication in QuarkDISFromSpline::TotalCrossSection(InteractionRecord) is preserved verbatim through this merge.
…ic modes Revert the dimuon-analysis-specific force-muonic convention (D+/D-/D0/D0Bar mu BR forced to 1.0, e and hadronic BRs forced to 0.0) introduced upstream 2025-04-09 and preserved through pavelzhelnin#3. Main MiaochenJin/SIREN ships PDG-physical BRs so downstream multi-cascade analyses can use SIREN charm without per-event reweighting by the inverse forced BR. Adds the Ds+ -> Hadrons e+ nu_e and Ds- -> Hadrons e- nu_ebar signature registrations omitted by Pavel (symmetric with the mu channels). The force-muonic mode remains available on pavelzhelnin/SIREN as a local patch for the dimuon analysis. Physical PDG values applied (CharmMesonDecay::BranchingRatio): - D+ / D-: e 0.1607, mu 0.176, hadrons 0.6633 - D0 / D0Bar: e 0.0649, mu 0.067, hadrons 0.8681 - Ds+ / Ds-: e 0.0654, mu 0.0654, hadrons 0.8692
# Conflicts: # projects/dataclasses/private/InteractionRecord.cxx # projects/interactions/CMakeLists.txt # projects/interactions/private/pybindings/interactions.cxx # projects/utilities/public/SIREN/utilities/Constants.h # python/SIREN_Controller.py
nickkamp1
left a comment
There was a problem hiding this comment.
This looks pretty good to go, I left a handful of inline comments, and I think we still need to address some of @austinschneider's comments. The biggest comments I have are:
- let's remove extraneous comments everywhere they appear. I saw a lot in this PR
- how do we want to handle building the pythia libraries for the PythiaDISCrossSection class? I think the code now assumes SIREN is built against Pythia libraries, but that is using a hard-coded path to Andy's pythia installation. we need to treat this more carefully for the main branch
There was a problem hiding this comment.
I think we do not want to ignore vendor updates, no?
| X(NuTau, 16) | ||
| X(NuTauBar, -16) | ||
| X(Charm, 4) | ||
| X(CharmBar, -4) |
There was a problem hiding this comment.
can we not use the c and cBar particles defined above? also probably good to group K0/K0Bar with the other mesons
There was a problem hiding this comment.
Done — removed X(Charm, 4) and X(CharmBar, -4); they were PDG-code duplicates of the X(c, 4) / X(cBar, -4)
entries defined above. Also moved X(K0, 311) / X(K0Bar, -311) next to X(K0_Short, 310) in the mesons section. The
only Charm/CharmBar callers were isQuark() (in Particle.cxx) and a switch case in
QuarkDISFromSpline::getHadronMass, both removed as dead code now that QuarkDIS emits D mesons directly rather than
bare charm quarks. The c / cBar entries (used by HNL and electroweak decay code via Constants::charmMass) stay as
the canonical charm-quark tag.
| double gen_prob = secondary_process_weighter_maps[idx].at(datum->record.signature.primary_type)->GenerationProbability(*datum); | ||
| physical_probability *= phys_prob; | ||
| generation_probability *= gen_prob; | ||
| // if (phys_prob == 0) { |
| } | ||
| } | ||
| inv_weight += generation_probability / physical_probability; | ||
|
|
| double target_prob = target_density * cross_section->TotalCrossSection(fake_record); | ||
| double total_xs = cross_section->TotalCrossSection(fake_record); | ||
| double target_prob = target_density * total_xs; | ||
| // if (total_xs == 0) { |
| if(signature == record.signature) { | ||
| // selected_prob += target_prob; | ||
| selected_final_state += target_prob * cross_section->FinalStateProbability(record); | ||
| double final_prob = cross_section->FinalStateProbability(record); |
There was a problem hiding this comment.
I think we want to add an optional install of Pythia to the SIREN build system, similar to how we have Marley handled now
There was a problem hiding this comment.
Agreed. I don't think we can realistically ship binary python wheels that include a full Pythia installation.
Gating this in CMake (as Nick is suggesting) is the best option for the moment. It will have to be only available for users willing to build from source themselves.
Side note:
I started to explore a concept where we would have a MARLEY (or other dependency) come in via external python package, but this can't be done purely via dynamic shared object loading at runtime on macOS and the other solutions are not very nice...
There was a problem hiding this comment.
That said, we should really consider making Pythia and MARLEY builds part of the CI testing.
There was a problem hiding this comment.
this logic probably needs to be checked against existing cross section classes... I don't think they all have TotalCrossSectionAllFinalStates. @MiaochenJin I thought we had agreed to update the way QuarkDISFromSpline handles TotalCrossSection(fake_record), no? i.e., it would check the final state and not double-count
There was a problem hiding this comment.
This was incorrectly pushed. Reverted now. The QuarkDISFromSpline class correctly handles rate calculation.
There was a problem hiding this comment.
This file's diff is now only whitespace changes, I think we can remove these changes.
…ontroller without changing constructor (PR125), please double check
…e end of text files
Co-authored-by: Austin Schneider <hogenshpogen@gmail.com>
Co-authored-by: Austin Schneider <hogenshpogen@gmail.com>
nickkamp1
left a comment
There was a problem hiding this comment.
One small comment in Constants.h, otherwise many comments in SIREN_Controller.py
| static const double DsMinusMass = 1.96835; // GeV | ||
| static const double BPlusMass = 5.27932; // GeV | ||
| static const double BMinusMass = 5.27932; // GeV | ||
| static const double CharmMass = 1.27; // GeV (charm-quark mass, meson-context legacy name) |
There was a problem hiding this comment.
@MiaochenJin can we not just use charmMass below? I don't think we should have both
There was a problem hiding this comment.
These are pretty extensive changes to SIREN_Controller. Are they all necessary? Inline comments below
| """ | ||
|
|
||
| # Define the primary injection process primary type | ||
| # Define the primary injection process primary type. |
There was a problem hiding this comment.
This change was not necessary before, and doesn't exist in the main branch, so I don't understand the comment about upstream PR 71. What has changed in this branch?
| sec_idist_list = list(secondary_injection_distributions[i_sec]) | ||
|
|
||
| # Add the position distribution | ||
| if fid_vol_secondary and self.fid_vol is not None: |
There was a problem hiding this comment.
why have we removed the fid_vol_secondary flag? This should be optional for the user
| :param list<dict<str,PhysicalDistribution> secondary_physical_distributions: List of dict of physical distributions for each secondary process | ||
| """ | ||
|
|
||
| # Define the primary physical process primary type |
| secondary_injection_process.primary_type = secondary_type | ||
|
|
||
| # Add the secondary position distribution | ||
| if fid_vol_secondary and self.fid_vol is not None: |
There was a problem hiding this comment.
same question, why are we removing this flag?
| assert(self.primary_injection_process.primary_type is not None) | ||
| # Use controller injection objects | ||
| self.injectors.append( | ||
| _injection.Injector( |
There was a problem hiding this comment.
where did these underscores come from? I think the pybind11 class is still Injector and not _Injector, and Weighter, not _Weighter
There was a problem hiding this comment.
These are monkey patched in the top-level SIREN __init__.py so that the interface is a little bit more convenient on the python side.
https://github.com/Harvard-Neutrino/SIREN/blob/main/python/__init__.py#L33-L45
…arm mass definition in constants.h
nickkamp1
left a comment
There was a problem hiding this comment.
all of my comments addressed with latest commit
Summary
Adds charm-DIS injection support to SIREN — primary neutrino-on-nucleon DIS
that emits D mesons directly as secondaries, plus the secondary D-meson
decay and energy-loss machinery needed to track D mesons through to their
final-state products. Used in IceCube high-energy charm analyses.
What this branch adds vs.
mainNew cross-section / interaction classes
QuarkDISFromSpline(projects/interactions/{public,private}/...QuarkDISFromSpline.{h,cxx}) —spline-based DIS cross section parameterised on charm-quark final states.
Loads upstream FITS splines (
differential+total), supportsCC / NC / GR
interaction_type_, and emits a[lepton, hadrons, D-meson]signature per primary × target × D-meson type. The D-meson set is selected
per-primary by parity (ν → {D⁰, D⁺}, ν̄ → {D̄⁰, D⁻}), guaranteeing charm
conservation in mixed-flavor instances.
PythiaDISCrossSection— Pythia8/LHAPDF-driven DIS cross section,same physics coverage as
QuarkDISFromSplinebut with Pythia-basedsampling. Drop-in alternative.
CharmMesonDecay+CharmMesonDecay3Body— semileptonic D-mesondecay (
D → K ℓ ν) and 3-body decay variants. Hadronization-classclasses (not Interaction); used as secondary collections.
DMesonELoss— D-meson energy-loss cross section, applies to Dmesons propagating through detector material as secondaries.
New resources / examples
resources/examples/example1/DIS_IceCube_charm.py— end-to-end charminjection example using SIREN's current Injector/Weighter idiom (no
SIREN_Controller, noAdd*Distribution(...)API). Runs to 100+ eventscleanly.
python/SIREN_Controller.py— back-compat patch for upstream'spython_interfacerefactor (PR New Python interface #71): updates to use_util.load_detector(experiment)and the new.distributionslistproperty in place of the removed
Add*Distributionsetters.New constants
projects::utilities::Constants— addsD0Mass = 1.86484(PDG D⁰),DPlusMass = 1.86966(PDG D⁺), and supporting decay-form-factorconstants (
alpha,m_star_squared, etc.).Build / dependency integration
projects/interactions/CMakeLists.txt— Pythia8 + LHAPDF detection andlink wiring, gated on user-provided
PYTHIA8_DIR/LHAPDF_DIRcachevariables (defaults configurable to a system install).
cmake/photospline_patches/+cmake/Packages/photospline.cmake—workaround for an upstream-pinned
photosplinerevision that's missingbsplvb_simple/bspline_deriv_nonzeroimpls. Patches restore the 3function bodies before
add_subdirectory(vendor/photospline). Whenupstream photospline is fixed, this dir can be deleted.
Tests
projects/interactions/private/test/CharmMesonDecay_TEST.cxx— unittests for D-meson decay kinematics (form factor, Q² distribution, decay
width normalisation, invariant-mass closure).
Verification status
cmake .. && make -j8) and Python pip install(
pip3 install -e .) both clean against Python 3.7 (CVMFS) and Python3.11 (Rocky 8) toolchains.
import siren; from siren import interactionssucceeds.DIS_IceCube_charm.py: exit 0,events written to disk.
QuarkDISFromSplineinstance withmixed-parity
primary_types = [NuMu, NuMuBar]returns{D⁰, D⁺}signatures forNuMuand{D̄⁰, D⁻}signatures forNuMuBar.tests/python/pytest suite (upstream's): 89 passed, 0 failed,exit 0.
Known limitations (intentionally out of scope, will land in follow-up PRs)
getIndices()leaveslepton_iduninitialised wheninteraction_type_ == 3(Glashow Resonance + charm). The path is notused in our physics analyses; cheapest fix is a
throwuntil the pathis properly supported.
CharmMesonDecayunimplemented decay modes print tostd::coutandsilently fall through. Should be
throws — behavioural change held backfor a separate PR.
cc̄pair is not modelled. The existingsplines and signatures encode the leading-charm (per primary parity)
convention only.
siren::interactions::CharmDISFromSpline+CharmHadronization:these classes existed in the pre-rebase form of this branch as an
intermediate "Charm quark → Hadronization → D meson" pathway. They are
fully superseded by
QuarkDISFromSpline's direct D-meson emission andhave been removed in this PR (they were dead code that never fired).