Prompted by @benjaminrose’s suggestions in this comment on the JOSS review:

Tests don't install the package via pip. They currently test the source code, good. But don't test what will be sent installed from pypi, less good. This explains how you are testing your project directory and not how it will be installed by users, https://hynek.me/articles/testing-packaging/. This is a valid test suite, but may not be doing what you expect.

Some background:
Since I didn’t have experience with Python packaging before working on sntools, I tried to stick close to the official recommendations (1, 2), which did not mention the src folder structure recommended in the linked blog post. (Though now that I’ve looked into it, it seems to me that this is hotly debated.)

I think my first step to improve packaging should be to set up a GitHub Action that automatically packages sntools and publishes it to PyPI whenever there’s a new release. (That will largely eliminate human error, which is probably the main source of packaging issues.) I’ll check how it’s possible to install & test as part of that action and will update this issue accordingly.

[based on discussions with Yixuan Jiang]

Some channels (most importantly es, where the contribution at low energies is strongest) currently enforce a Cherenkov threshold on the outgoing particle. While this makes sense for water Cherenkov detectors (to avoid simulating undetectable events), we should remove this threshold for LS detectors.

Following the end of support for Python 2.7 and in line with many other packages in the scientific Python ecosystem (including sntools’ dependencies, NumPy and SciPy), I plan to drop support for Python 2.7 soon.

• v0.7.4 (to be released soon) will print a warning message when running under Python 2.7 and point people here.
• If any further bugfix updates (v0.7.x) appear, they will remain compatible with Python 2.7
• Future feature updates (either v0.8 or v1.0) will drop support for Python 2.7 and will likely require Python 3.6 or 3.7.

This issue will be updated with any progress on this.

These interactions could make up about 5% of events for a SN burst detection in Super-K. A theoretical calculation of the cross-sections was performed by Langanke et al., arXiv:nucl-th/9511032; experimental data appears to be non-existent.

The main issues here are

• the low energy of events (typically gammas with just above 5 MeV)
• large uncertainty of the cross-section

and as a result, this interaction is not currently considered in the HK Design Report.

snewpy is a SNOwGLoBES frontend developed by members of the SNEWS 2.0 collaboration. Like SNOwGLoBES, it’s not a replacement for sntools; but it does have features that would be useful in sntools, too:

• support for more flavour transformations
• support for more supernova models

Time- and energy-independent flavour transformations should be straightforward to add to sntools, analogous to the current implementation of MSW effect in normal/inverted mass ordering. (However, it will be necessary to change the command line interface to e.g. --transformation=AdiabaticMSW_NMO (instead of --ordering=normal), once additional transformations are available.)

Time- or energy-dependent flavour transformations will require further changes to sntools.

For additional supernova models, I need to check if it’s possible to write an adapter that maps snewpy’s interface to the one expected of an input format in sntools.

Dear Jost,

I've noticed a minor typo in the sntools/format/init.py file that affects compatibility with the latest version of snewpy. Specifically, at line 165, the correct code should be:

sn_model = getattr(import_module('snewpy.models.ccsn_loaders'), format)(abspath(file))

Could you modify?

Best regards,
Pierre-Alexandre

Start time and duration are currently hard-coded in ibd.py to go from 15ms to (15+335=)350ms. https://github.com/JostMigenda/sntools/blob/01d3ded3afa2f0c2a991d55f4e42e36f1491a5ab/ibd.py#L223
https://github.com/JostMigenda/sntools/blob/01d3ded3afa2f0c2a991d55f4e42e36f1491a5ab/ibd.py#L228

Instead, we could read the start & end time from the input file.

At the proto-collaboration meeting, Neil McCauley pointed out that nu_e CC interactions on oxygen-17/18 have a much lower threshold energy than the analogous oxygen-16 interaction, which leads to much higher cross-sections at low energies.

After a brief search, I’ve found two relevant papers:

• Haxton: PRD 36 (1998) 2283 calculates the cross-sections theoretically. Fig. 1 shows the cross-section; the "bump" in the solid line at low energies comes from oxygen-17/18. This implies that those isotopes are negligible (compared to oxygen-16) at energies above ~20 MeV.
• Haxton & Robertson, arXiv:nucl-th/9806081 calculates the cross-sections averaged over the solar boron-8 flux and finds that the oxygen-18 event rate is about 6.7% of the total solar neutrino signal (i.e. ~14× smaller than elastic scattering on electrons), while the oxygen-17 interaction is negligible (cross-section and natural abundance are about 3.5 and 5 times smaller than for oxygen-18).

My very rough guess (based on this and the HK Design Report) is that we might see 100–200 events from this oxygen-17/18 channel in HK. This would be a little more than the oxygen-16 channel interaction in the worst case (Totani1998 model, no oscillations), but an order of magnitude fewer events then we’d see for the oxygen-16 channel in the NH/IH scenarios.
Also, the oxygen-17/18 channel would depend only weakly on the mean neutrino energy.

The NuHepMC output format (see arXiv:2310.13211) that’s been in the works for a while now looks like it’s almost ready. Notably, the new Hyper-K detector simulation expects to use this as the standard input format for all event generators; so sntools should add support for this output format.

See also the NuHepMC organisation on GitHub.

Various minor documentation improvements, as suggested by @benjaminrose in this post in the JOSS review.

• Make clear(er) that sntools is a command line app.
• Compiled version of documentation.
• First example should run more quickly.
• Simplify arguments of first example

In addition to being an event generator, sntools can be used in other code as a library that implements various cross sections. (See this SNOwGLoBES PR for an example.)
This should be added to the documentation.

This paper by Ricciardi et al. (2022) has an improved evaluation of the IBD cross section, showing up to ~1% (at 100 MeV) difference compared to the Strumia & Vissani XS currently implemented, as well as a more precise evaluation of the uncertainty.

It would be nice to implement this updated cross section.

Currently, one can either generate events for all supported channels in the respective detector (default) or for a single channel (with e.g. the --channel ibd command line argument). It may be useful to select multiple interaction channels, e.g. by allowing multiple occurrences of the argument, e.g. --channel o16e --channel o16eb.

See eqn. (14) in the Strumia/Vissani paper.

For reproducibility, it'd be nice to have the option to set a seed in genevts.py

From our discussions, you've noted that you used both random and np.random So these should both have seeds sets. Whether they're the same, indpendently settable, or have an offset (e.g. SEED and SEED + 10) I have no opinion on

These interaction channels are mainly relevant at low energies due to their energy threshold, which is much lower than that of equivalent interactions on carbon-12.

However, carbon-13 makes up only about 1% of carbon nuclei in nature. Thus, even at low energies I’d estimate that the event rate of these interactions is about an order of magnitude smaller than that of neutrino-electron scattering. That would make this a very minor correction.

A theoretical calculation (including a nice numerical fit) of the cross-section is available.

Discovered this issue while working on a PR for SNOwGLoBES: The oxygen-16 CC cross sections are implemented using a fit from arXiv:1809.08398, that has four "excitation energy" parameters. If the neutrino energy is barely higher than one of these excitation energies, the cross section has an unphysical spike (because the fit function diverges as the difference between neutrino energy and excitation energy goes to zero).

Sample code to demonstrate:

from sntools.interaction_channels import o16e
for i in range(1,5):
    eG, a, b, c = o16e.fit_parameters[i]
    eNu = 29.350001 # excitation energy of the 4th group: 29.35 MeV
    print(f"*** {i}: eNu={eNu}, eG={eG} ***")
    if eNu - eG > 0:
        d = log10(eNu**0.25 - eG**0.25)
        log_sigma = a + b * d + c * d**2
        print(f"d = {d}\nlog_sigma = {log_sigma}\npartial_dSigma_dE(eNu, eNu-eG, {i}) = {o16e.partial_dSigma_dE(eNu, eNu - eG,i)}\n")

In reality, the difference between eNu and eG needs to be above the detector threshold (e.g. 0.8 MeV in water Cherenkov detectors). Enforcing that may be sufficient to fix this?